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A  feedback  controller  design  method  has  been  formulated  which  applies  to  linear 
time  varying  systems.  The  motivation  behind  this  technique  is  the  design  of  an 
autopilot  for  an  air-to-air  missile.  The  missile  under  study  is  the  extended  medium 
'  range  air-to-air  technology  (EMRAAT)  missile,  a  theoretical  bank-to-turn  missile 
which  is  under  study  by  the  United  States  Air  Force.  Knowledge  gained  by  this 
missile  will  be  applied  to  future  bank-to-turn  missiles. 

Conventional  autopilot  design  techniques  use  pole  placement,  linear  quadratic 
regulators,  linear  quadratic  Gaussian/loop  transfer  recovery  techniques,  or  eigen- 
structure  assignment  methods.  These  techniques  are  applied  to  a  linear  model  which 
depends  on  time  varying  flight  parameters.  In  applying  these  methods,  the  false  as- 
sumption is  made  that  the  linear  model  is  slowly  varying.  Since  the  system  is  changing 
rapidly,  no  theoretical  basis  exists  to  support  the  success  of  the  resulting  controller. 
The  proposed  design  method  finds  state  feedback  gains  which  cause  the  closed  loop 
linearized  system  to  be  stable  with  respect  to  a  specified  Lyapunov  function.  Even 
though  flight  parameters  change  rapidly,  local  stability  around  the  operating  point  is 

achieved. 

iv 


The  design  algorithm  finds  feedback  gains  which  place  the  eigenvalues  of  the 
derivative  of  a  given  Lyapunov  function.  Limitations  on  placement  of  these  eigenval- 
ues are  stated.  This  concept  is  applied  to  a  second  order  linear  time  varying  system 
where  ordinary  pole  placement  techniques  fail.  The  design  method  is  then  applied 
to  a  linearized  model  of  the  EMRAAT  missile  which  is  a  function  of  time  varying 
flight  parameters.  Feedback  gains  are  generated  as  a  function  of  these  flight  param- 
eters. It  is  discovered  that  the  gains  depend  on  dynamic  pressure,  mach  number, 
and  angle- of- at  tack.  A  combination  of  interpolation  and  polynomial  fitting  is  used 
to  create  a  look-up  table  for  the  gains.  The  resulting  controller  is  programmed  into 
a  nonlinear  simulation  which  runs  missile  and  target  scenarios.  Small  miss  distances 
are  achieved. 


CHAPTER  1 
INTRODUCTION 

The  original  goal  of  this  dissertation  was  to  formulate  an  autopilot  design  method 
which  stabilizes  the  flight  of  a  bank-to-turn  air-to-air  missile.  The  result  was  the  de- 
velopment of  a  controller  design  method  which  applies  to  time  varying  linear  systems. 
When  used  on  a  nonlinear  system,  local  stability  is  achieved.  The  remainder  of  this 
chapter  gives  a  brief  description  of  the  missile  used  in  this  study,  a  summary  of  related 
works  preceding  this  study,  and  the  purpose  and  outline  of  this  dissertation. 

An  air-to-air  missile  is  launched  from  a  military  aircraft  with  the  intent  of  in- 
tercepting an  enemy  aircraft  or  target.  Bank-to-turn  (BTT)  missiles  have  wings 
giving  them  greater  maneuverability  than  the  conventional  skid-to-turn  (STT)  mis- 
siles which  have  fins.  Thus,  the  target  must  work  much  harder  to  evade  a  BTT 
missile.  The  dynamics  of  BTT  airframes  are  highly  unstable,  making  for  a  difficult 
controls  problem. 

The  missile  under  study  in  this  paper  is  the  extended  medium  range  air-to-air 
technology  (EMRAAT)  missile  as  shown  in  Figure  1.1  [l].  The  EMRAAT  missile 
is  a  theoretical  missile  formulated  by  the  Air  Force  with  the  intent  of  studying  the 
feasibility  of  the  BTT  concept.  The  knowledge  gained  from  the  study  of  this  missile 
will  be  used  by  the  Air  Force  if  it  ever  decides  to  make  a  real  BTT  missile. 

The  EMRAAT  missile  is  equipped  with  a  seeker.  If  the  seeker  is  infrared  based, 
then  it  measures  the  line  of  sight  angles  to  the  target.  In  the  case  that  the  seeker  is 
radar  based,  then,  in  addition  to  the  line  of  sight  angles,  the  range  and  range  rate 
of  the  target  are  measured.  This  information  is  passed  to  the  guidance  law  which 
determines  the  desired  normal  accelerations  needed  to  intercept  the  target.     The 
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Figure  1.1.  The  EMRAAT  missile  [1]. 


autopilot  attempts  to  achieve  these  accelerations  by  applying  the  proper  control  to 
the  control  surfaces. 

1.1     Earlier  Works 

Early  autopilot  designs  separated  the  missile  dynamics  into  two  or  three  lower 
order  models.  These  models  were  linearized,  and  classical  techniques  were  used  to 
stabilize  them.  Other  designs  applied  more  advanced  techniques  to  higher  order  lin- 
earized models.  Some  of  these  techniques  are  pole  placement,  linear  quadratic  regula- 
tors, linear  quadratic  Gaussian/loop  transfer  recovery  techniques  and  eigenstructure 
assignment  methods  [2]  Feedback  controllers  would  be  designed  as  a  function  of  the 
linearized  models.  Gains  would  then  be  scheduled  against  the  flight  parameters  on 
which  the  models  depend.  Gain  scheduling  is  a  popular  technique  used  in  control 
designs  and  has  been  studied  in  depth  by  Shamma,  Athans,  and  Cloutier  [3,4,5,6]. 

The  problem  with  these  design  methods  is  that  they  assume  that  the  system  is 
time  invariant  or  slowly  varying.  The  linearized  models  depend  on  rapidly  changing 
flight  parameters.  Although  these  methods  often  work,  no  theoretical  basis  exists 
to  support  the  success  of  these  designs.  More  recently,  H^  and  jj,  synthesis  design 
techniques  have  been  used  to  formulate  one  dynamic  controller  which  would  stabilize 
the  system  for  all  modeled  uncertainties  [7,8].  A  disadvantage  of  H^  controllers  is 
that  they  have  very  high  orders.  Also,  the  existence  of  a  robust  H^  controller  which 
satisfies  the  performance  requirements  is  not  guaranteed.  The  only  way  to  determine 
the  stability  of  a  given  design  is  to  test  it  using  a  nonlinear  simulation. 

Desoer  [9]  states  if  a  given  time  varying  system 

X  =  A(t)x  (1.1) 

is  stable  for  each  t  when  t  is  frozen,  then  it  is  not  necessarily  stable  when  t  is  allowed 
to  change.    However,  (1.1)  is  stable  if  A(t)  changes  slowly.    A  conservative  upper 
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Figure  1.2.  Angle  of  attack  from  one  trajectory  of  the  EMRAAT  missile. 

bound  on  supt>0  \\A(t)\\  is  given  which,  if  enforced,  guarantees  asymptotic  stability 
of  (1.1).  As  will  now  be  shown,  this  limit  is  too  restrictive  for  the  EMRAAT  missile. 
Figure  1.2  shows  the  angle  of  attack  of  the  trajectory  of  a  missile  flown  using  an 
already  existing  autopilot.  Desoer  proposes  the  Lyapunov  function 


V(x,i)  =  xT(ei/  +  P(i))x 


'1.21 


where  e%  >  0  and  P(t)  is  chosen  so  that 


(T 


A1{t)P{t)  +  P{t)A{t)  =  -ZI 


1.3) 


Desoer  gives  a  bound  on  V, 


V  <  xTx[-3  +  2£laM  + 


where  the  following  definitions  are  made: 


(1.4) 


cim  '■=  sup  ||/4.(i)||  <  oo 


'1.5] 


(To  is  positive  and 

Re  Xi[A(t)\  <  -2(7o  Vi,  V<  >  0  (1.6) 

m  is  a  constant  and  depends  only  on  a0  and  aM  so  that 

||exp(rA(i))||<mexp(-<70r),  Vr  >  0  Vi  >  0  (1.7) 

Also, 

ajfcf  :=  sup  ||A(t)||.  (1.8) 

i>0 

If  e1  is  allowed  to  be  very  small,  then  from  (1.4),  stability  would  result  if 

4<r2 
°m  <  —^  (1.9) 

For  the  linear  pitch  model  of  the  given  trajectory,  a0  =  15.75.   At  .04  seconds  into 
the  trajectory  the  closed  loop  matrix  is 


A  = 


-3.2533       .8997 
-1678.2    -83.50 


(1.10) 


At  this  instant  in  time,  m  =  22.88  and  hM  =  2231.  jjj  =  .0036.  From  this  we  see 
that  the  linear  system  is  changing  much  faster  than  the  limit  shown  in  (1.9).  The 
system  is  not  slowly  varying.  Note  that  m  was  only  found  for  t  =  .04  sec.  and  not 
for  all  time.  If  a  greater  m  were  found  for  the  rest  of  the  trajectory,  then  the  limit  in 
1.9  would  be  even  more  restrictive  and  the  result  would  be  the  same.  This  test  does 
not  show  that  the  system  is  unstable.  In  fact,  a  Lyapunov  function  is  known  which 
shows  that  this  linear  system  is  stable  but  Desoer's  theory  can  not  support  this  fact. 
Wilson,  Cloutier,  and  Yedavalli  [1]  give  a  stability  condition  for  a  constant  linear 
system  with  time  varying  uncertainties.  Given  a  system, 

X=[A  +  E(t)]x  (l.H) 

if  a  positive  definite  P  can  be  found  which  solves  the  Lyapunov  equation, 

PA  +  ATP  +  2I  =  0  (1.12) 
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then  the  system  will  be  asymptotically  stable  if 
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where  amax  is  the  maximum  singular  value.  This  idea  is  initially  interesting  because 
it  requires  the  error,  E(t),  to  be  bounded,  but  does  not  constrain  the  time  variation. 
The  above  result  can  not  be  easily  applied  to  the  EMRAAT  missile.  This  will  be 
shown  by  applying  the  result  to  a  trajectory  of  the  EMRAAT  missile  whose  angle  of 
attack  is  shown  in  figure  1.3.  In  order  to  check  the  stability  of  the  EMRAAT  missile. 
A  must  be  chosen  since  it  may  be  any  matrix  which  is  stable.  It  is  believed  that  the 
best  choice  of  A  should  come  from  a  closed  loop  matrix  in  the  trajectory.  With  this 
in  mind,  A  was  taken  at  t  =  2  sec  into  the  trajectory. 


A  -  BK 


-4.477       .9526 
-1.684    -79.47 


1.14) 
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Figure  1.4.  A  plot  of  <rmax[E(t)]  versus  time. 

E(t)  is  the  result  of  subtracting  this  matrix  from  the  closed  loop  system  from  the 
rest  of  the  trajectory.  Solving  for  P  in  (1.12)  gives 


P  = 


17.29 
-.0454 


-.0454 
.012 


(1.15) 


Computing  the  upper  bound  in  (1.13)  gives 

1 


a. 


,(P) 


=  .0578 


(1.16) 


Figure  1.4  shows  the  plot  of  <rmax[E{t)}.  Obviously  since  <rmax[E{t))  is  greater  than 
.0587  for  most  of  the  trajectory  then  (1.13)  is  not  satisfied.  This  procedure  was 
repeated  by  taking  the  constant  closed  loop  matrix  from  every  point  in  the  trajectory, 
and  the  result  was  the  same.  The  time  varying  linear  system  from  this  trajectory  is 
known  to  be  stable  because  a  Lyapunov  function  exists  which  can  show  this.  The 
upper  bound  in  (1.13)  can  not  be  easily  used  if  at  all  to  support  the  claim  of  stability. 
Desoer  [9]  gives  a  stability  limit  on  the  rate  of  variation  of  a  given  closed  loop 
time  varying  system.    Wilson,  Cloutier,  and  Yedavalli  [1]  give  a  stability  limit  on 
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the  size  of  a  time  varying  uncertainty.  Both  of  these  tests  were  found  to  be  too 
conservative  when  applied  to  an  already  existing  stable  autopilot,  and  the  results 
were  inconclusive.  In  addition,  these  tests  are  analytic  tools.  As  yet,  no  design 
procedure  exists  which  can  give  a  closed  loop  linear  system  that  will  pass  either 
test.  It  would  be  desirable  to  formulate  a  design  method  that  can  yield  a  closed 
loop  time  varying  linear  system  which  satisfies  a  less  conservative  stability  condition. 
Vidyasagar  [10]  gives  two  important  results.  The  first  can  also  be  found  in  Hahn  [11] 
and  says  that  given  a  time  varying  system 

X  =  A(t)x  (1.17) 

if  a  positive  definite  matrix  P  can  be  found  so  that  the  matrix 

PA(t)  +  AT(t)P  (1.18) 

is  negative  definite  for  all  time,  then  the  system  is  asymptotically  stable. 

The  second  result  allows  the  first  result  to  be  applied  to  a  nonlinear  system.  Given 
a  nonlinear  system  of  the  form 

x  =  /(t,x) 
where 

/(i,0)  =  0 

and  /  is  continuously  differentiate,  then  let 

df(t,x) 


A(t)  = 


x=0 


dx 
and  assume  that 

limsu  „  !!/(«,*) -A(Qx|| 


sup 
WHO  t>5 


and  A  is  bounded.  If  0  is  an  asymptotically  stable  equilibrium  point  of  i  =  A(t)z  for 
all  time,  then  0  is  a  locally  stable  equilibrium  point  for  the  system 

x  =  /(t,x) 

These  ideas  can  be  applied  in  the  following  way.  Given  a  nonlinear  system 

x  =  /(x,w,u)  (1.19) 


define 

./.._*       df  df 

(1.20) 


A(x,w)=| 


S(x,w)-|£ 

^jkj  »v  Jnominal 


\X->™  )  nominal 

where  x  is  a  vector  containing  the  states  and  w  is  a  vector  containing  additional 
system  parameters.  At  regions  near  the  operating  point,  the  system  becomes  close 
to 

Ax  =  A(x,w)Ax  +  5(x,w)Au  (1.21) 

where  Ax  and  Au  are  small  perturbations  between  the  states  and  inputs  and  the 
operating  point.  We  would  like  to  find  a  feedback  control  law, 

u  =  -iT(x,w)x  (1-22) 

so  that 

P[A(x,  w)  -  B(x,  w)#(x,  w)]  +  [A(x,  w)  -  B(x,  w)A'(x,  w)}TP  (1.23) 

is  negative  definite  for  all  x  and  w  where  P  is  positive  definite.  If  such  a  control  law 
is  found,  then  the  perturbations  from  the  nominal  trajectory  are  locally  stable  for 
the  system 

x  =  /(x,w,  -A'(x,w)x)  (1.24) 

Shahruz  and  Behtash  [12]  give  one  control  law  which  places  some  of  the  eigenval 
of  (1.23)  for  the  case  where  P  =  I.   However  unnecessary  limitations  on  where  th 
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eigenvalues  can  be  placed  exist.  Also,  while  many  feedback  control  laws  make  (1.23) 
negative  definite,  Shahruz  and  Behtash  give  only  a  small  subset  of  these  control  laws. 

1.2     Purpose 

This  paper  gives  an  algorithm  which  can  stabilize  a  linear  time  varying  system  by 
placing  the  eigenvalues  of  (1.23)  between  desired  bounds  and  has  fewer  limitations 
than  the  control  law  in  Shahruz  and  Behtash  [12].  Limitations  on  the  placement  of 
these  bounds  are  stated.  This  algorithm  applies  to  all  systems  for  which  a  constant 
positive  definite  P  exists  so  that  (1.23)  can  be  made  negative  definite  for  all  time. 
First,  theory  is  presented  leading  to  the  formulation  of  this  algorithm.    Then,  to 
demonstrate  the  usefulness  of  this  algorithm,  it  is  applied  to  a  linear  time  varying 
system  where  normal  pole  placement  techniques  fail.  This  algorithm  is  then  applied 
to  the  EMRAAT  missile.   A  nonlinear  model  is  made,  and  from  this,  a  time  vary- 
ing linearized  model  is  generated  which  is  a  function  of  several  flight  parameters. 
Gains  are  computed,  and  their  dependence  on  flight  parameters  determined.  A  gain 
scheduling  scheme  is  then  implemented  yielding  local  stability  around  the  operating 
point.  The  resulting  design  is  tested  in  a  nonlinear  simulation.  Finally,  the  results  of 
this  test  are  given,  and  the  usefulness  of  this  design  technique  is  evaluated. 


CHAPTER  2 
THEORY  BEHIND  THE  DESIGN  METHOD 

The  material  in  this  section  gives  the  theory  leading  to  the  proposed  controller 
design  method.  The  first  section  presents  a  geometric  interpretation  of  Lyapunov's 
linear  stability  theorem.  The  effects  of  control  on  the  velocity  field  of  a  system  is 
studied  in  the  next  section.  The  nonemptiness  and  convexity  of  the  set  of  all  feedback 
gams  which  stabilize  a  system  with  respect  to  a  given  Lyapunov  function  are  then 
discussed.  Finally,  an  iterative  procedure  that  finds  one  element  of  this  set  is  given. 

2T — A  Geometric  Interpretation  of  Lvapunov's  Linear  Stability  Theorem 

To  understand  why  time  invariance  is  a  necessary  assumption  for  eigenstructure 
design  methods  the  following  second  order  linear  time  varying  system  was  studied. 
The  example  given  now  is  from  Vidyasagar's  example  5.3,109  [10]  and  can  also  be 
found  in  Khalil  [13].  Given  the  following  system 


A(t)x 


where 


£U)  =  -I  +acos2(t)     1  -  asin(t)cos(t) 

-1  -  a  sin(i)  cos(i)  -1  +  a  s'm2(t) 

Vidyasagar  [10]  notes  that  the  transition  matrix  is  given  by 


*(*,0). 

and  the  characteristic  equation  is 


o-l  < 


cos(t)     e  fsin(i!) 
-e^-^sinOO    e-*cos(<) 


(2.i; 


(2.2) 


(2.3) 


A2  +  (2-a)A  +  (2-a)  =  0 


(2.4) 


The  roots  of  (2.4)  have  negative  real  parts  for  1  <  a  <  2.  The  exponents  in  the  first 
column  of  $(i,  0)  indicate,  however,  that  the  system  is  unstable  for  these  values  of  a. 
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Figure  2.1.  The  trajectory  of  (a)  the  system  (2.1)  and  (b)  a  frozen  system. 

This  system  was  simulated  using  MATLAB  with  a  =  1.5  and  x0  =  [l,0]r.  The 
eigenvalues  of  A(t)  are  -l  +  &j  and  -f-f  j.  Figure  2.1a  shows  the  state  trajectory 
of  the  system  2.1  where  Xl  is  assigned  to  the  horizontal  axis,  and  x2  is  assigned  to 
the  vertical  axis.  This  plot  demonstrates  the  instability  of  the  system.  If  the  system 
were  frozen,  i.e.  A(t)  =  A(0),  then  the  stable  trajectory  of  Figure  2.1b  would  result. 
This  trajectory  is  shown  with  the  velocity  vector  field  of  the  frozen  system. 

To  explain  the  instability  in  Figure  2.1a,  we  will  now  discuss  the  time  varying 
velocity  field  of  equation  (2.1).  The  plots  in  Figure  2.2  show  the  trajectory  of  the 
time  varying  system  in  the  state  plane  for  eight  instances  in  time.  For  each  plot,  the 
velocity  vector  field  for  that  instant  in  time  is  superimposed  on  to  the  trajectory.  Each 
plot  in  Figure  2.2  shows  the  boundaries  of  four  pie-shaped  regions  which  this  paper 
refers  to  as  positive  and  negative  regions.  Two  positive  regions  exist  that  are  defined 
as  the  set  of  all  points  whose  velocities  have  an  outer  radial  component,  i.e.  each  of 
these  velocities  have  a  component  pointing  directly  away  from  the  origin.  Likewise, 
there  are  two  negative  regions  containing  velocities  with  inner  radial  components. 
The  plots  show  that  the  boundaries  separating  the  positive  and  negative  regions 
rotate  in  a  clockwise  direction.  Also,  the  current  position  in  the  trajectory  remains 
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inside  one  of  the  positive  regions.  Since  the  velocity  of  the  system  always  has  an 
outer  radial  component,  the  magnitude  of  the  state  vector  is  always  increasing  so 
that  the  system  (2.1)  is  unstable. 

A  direct  relationship  exists  between  the  positive  and  negative  region  boundaries 
and  the  eigenvectors  of  the  system.  In  this  case,  the  eigenvectors  are  rotating  clock- 
wise at  the  same  angular  velocity  as  the  region  boundaries.  It  can  be  shown  that  if  a 
system  has  constant  eigenvectors  with  time  varying  eigenvalues  which  have  negative 
real  parts  for  all  time  then  the  system  is  stable.  The  stability  question,  however,  is 
not  as  easy  to  answer  for  the  case  with  moving  eigenvectors. 

If  the  system  has  no  positive  regions,  then  it  will  be  stable  regardless  of  how 
quickly  it  changes.  However  the  converse  is  not  true  in  general.  Figure  2.3  shows  a 
state  vector  and  its  corresponding  velocity  for  some  dynamic  system.  The  angle  9 
can  be  computed  in  the  following  way. 

T  • 

cos(d)  =  nT^T  (2.5) 


xx 


For  x  to  have  an  inner  radial  component  then 


2  ~      -    2 


This  is  true  when 

xTx<  0 
where 

x  =  Ax 


(2.6) 


(2.7) 


(2.8) 


So  if  x  Ax  is  negative  everywhere  for  all  time,  then  the  system  is  stable.  Since  the 
system  is  linear,  it  is  sufficient  to  check  the  sign  of  xTAx  on  the  unit  circle.  If  all 
velocities  on  the  unit  spheroid  point  inside  that  boundary,  then  the  same  is  true  for 
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Figure  2.2.  The  trajectory  and  velocity  field  of  (2.V 
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Figure  2.3.  x  and  x 

all  spheroids  centered  on  the  origin.   This  leads  to  the  following  stability  condition 
which  is  given  without  proof. 

Stability  Condition  1:  If  xTx  =  xTAx  <  0  V  x  :  xTx  =  1,  V  t  >  0  for 
the  system,  x  =  A(t)x,  then  the  origin  is  a  global  asymptotically  stable 
equilibrium  point. 

Figure  2.4  illustrates  an  example  of  a  second  order  system  which  meets  Stability 
Condition  1.  Note,  that  while  all  figures  given  so  far  represent  second  order  linear 
systems,  the  above  stability  condition  applies  to  linear  systems  of  any  order. 

The  above  stability  condition  can  be  made  less  conservative  by  expanding  the 
class  of  shapes  to  ellipsoids  defined  in  the  following  way. 


xTPx  =  1 


(2.9) 


where  P  is  a  symmetric,  positive  definite  matrix.  Figure  2.5  shows  an  ellipse  of  th 
form  (2.9)  defined  for  a  second  order  system.  The  state  vector,  x.  is  drawn  to 


some 
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Figure  2.4.  An  example  of  a  second  order  system  with  negative  radial  velocity 
ponents 


com- 


Figure  2.5.  One  velocity  vector  on  the  ellipse  xTPx  =  1 
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point  on  the  ellipse,  and  its  velocity,  x,  is  also  shown.  The  normal  to  the  ellipse  at  x 
is  Px.  To  ensure  stability,  it  is  sufficient  to  require  that  the  projection  of  x  onto  the 
normal  of  the  ellipse,  Px,  is  negative  everywhere.  The  resulting  term,  xTPx,  is  the 
normal  velocity  component  on  the  ellipse.  The  main  thrust  of  this  paper  is  to  be  able 
to  control  x  Px  and  to  make  this  normal  velocity  component  negative  everywhere 
for  all  time  in  order  to  achieve  stability.  As  before,  the  linearity  of  the  system  allows 
one  to  only  check  xTPx  on  the  unit  circle.  This  generalizes  the  previous  stability 
condition  to 

Stability  Condition  2:  For  the  system  x  =  A(t)x,  if  xTPx  =  xTPA(t)x  < 
0   Vx:x   x  =  1,    V  i  >  0  for  some  constant  positive  definite  matrix  P, 
then  the  origin  is  a  global  asymptotically  stable  equilibrium  point. 

The  above  condition  applies  to  linear  systems  of  any  order. 

Lyapunov  [14]  stated  that  if  any  positive  definite  function  of  the  system  states 
V(x)  is  always  decreasing,  i.e.  V(x)  <  0,  then  the  system  is  stable  (asymptotically 
stable  for  linear  systems)  [10].  For  the  linear  case,  let 

V(x)  =  xTPx  (2.10) 

V(x)  is  positive  definite  if  and  only  if  P  is  a  positive  definite  matrix.    Taking  the 
derivative  gives 

V(x)    =    xTPx  +  xTPx  (2.H) 

=    xTATPx  +  xTPAx  (2.12) 

=    xT[ATP  +  PA]x  (2.13) 

Lyapunov's  stability  condition  becomes 

xT[ATP  +  PA]x  <  0  V  x,  V  t  >  0  (2.14) 
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Si 


nice 


xTPx  =  xTPAx  =  -xT[ATP  +  PA]x  (2.15) 

Stability  Condition  2  is  equivalent  to  Lyapunov's  linear  stability  condition. 

It  would  be  useful  to  compute  the  lower  and  upper  bound  of  the  rate  of  decay  of 
a  given  positive  definite  function, 

V(x)  =  xTPx 
of  the  states  of  a  time  varying  linear  system  by  evaluating, 

[^rmn\[PA(t)  +  AT{t)P] 
[^maX-[PA(t)  +  AT(t)P] 


(2.16) 


c  :=  mm 

t 


c  :=  max 

t 


So, 


c  <  -xT[PA(t)  +  AT(t)P]x  <  c  V  x  :  xTx  =  1,    V  t 


.(2.17) 
(2.18) 

(2.19) 


Then  c  and  c  are  the  lower  and  upper  bounds  respectively  of  the  normal  velocity 
components,  xTPx,  on  the  unit  spheroid.  If  c  <  0,  then  the  system  is  stable. 

2J2 Adding  Control  to  the  System 

The  result  in  section  2.1  is  an  analytical  tool  only.  A  design  procedure  is  needed 
to  control  stability  for  the  system 


x  =  A(t)x  +  B(t)u 


(2.20) 


where  u  is  the  control  vector  and  x  is  the  state  vector.  This  section  focuses  on  two 
questions: 

1)  What  effect  does  u  have  on  the  velocity  field  of  (2.20)? 

2)  Can  the  normal  velocity  component  xTPx  on  the  spheroid  xTx  =  1 
be  controlled? 


19 


-i 

-IS 

-2 


•1J        -I         -0.5         0         0.5  1  IS 

xl 


7 

1 

■ 

B 

0.5 

x£s*£=Ax 

■a 

0 

\           — ' 

-05 

1 

-1 

fll 

-1-5 

I 

-1 

-2 

-1.5 

-OS 

0 

OS 

1          1J         2 

xl 

Figure  2.6.  A  velocity  vector  (a)  without  control  and  (b)  with  control. 

To  answer  question  1),  consider  a  second  order  single  input  system  frozen  at 
instant  in  time. 


some 


x  =  Ax  +  Bu 


(2.21] 


A  is  a  2  x  2  matrix,  and  B  is  a  two  dimensional  column  vector,  u  is  a  scalar  input 
and  can  take  on  any  real  value.  If  u  =  0  then  the  velocity  field  of  the  system  is 
x  =  Ax.  Figure  2.6a  shows  the  velocity,  x  =  Ax,  of  some  state  vector,  x,  in  the 
system.  The  control  vector,  5,  is  also  shown.  When  u  f  0,  x  has  the  extra  term, 
Bu.  The  direction  of  Bu  is  constant,  but  its  magnitude  is  directly  proportional  to  u. 
Figure  2.6b  shows  many  possible  values  for  x  by  sweeping  u  through  a  wide  range  of 
values  in  small  increments.  As  this  diagram  shows,  the  arrow  head  of  each  velocity 
vector  can  be  placed  on  the  line  drawn  parallel  to  B  and  intersecting  the  arrow  head 
of  Ax.  This  demonstrates  that  velocities  can  be  controlled  along  the  space  spanned 
by  the  columns  of  B.  The  following  theorem  answers  question  2). 
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Figure  2.7.  A  case  in  which  the  normal  velocity  component  cannot  be  controlled. 
Theorem  1 

Consider  the  system,  x  =  A(t)x  +  B(t)\l,  where  dim(x)  =  n,  dim(u)  =  m 
and  A(t)  and  B(t)  have  compatible  dimensions.  Let  P  be  a  constant 
positive  definite  matrix.  The  normal  velocity  component,  xTPx,  can  be 
arbitrarily  set  to  any  value  with  the  right  choice  of  u  ,  at  a  given  time  t, 
at  any  point  x,  on  the  spheroid,  xTx  =  1,  unless  x  is  contained  in  the  set 

S(t)  :=  p-lspan(B±{t))  n  {x  :  xTx  =  1}  (2.22) 

where  BL(t)  is  a  basis  of  column  vectors  orthogonal  to  span(B(t)).  If 
x  €  S(t)  then 

xTPx  =  xTPA(i)x  (2.23) 

and  this  velocity  component  cannot  be  controlled  at  time  t. 


Theorem  1  gives  the  parts  of  the  unit  spheroid  for  which  the  velocity  component 
x    Px  is  uncontrollable  and  can  therefore  be  used  to  determine  if  a  linear  time 
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system  is  stabilizable  with  an  appropriate  choice  of  constant  positive  definite  P.   If 
there  exists  a  constant  positive  definite  P  so  that 

mfX^)  xTpA^x  <  °  (2.24) 

the  system  is  stabilizable.  The  above  condition  is  equivalent  to  requiring  the  following 
matrix  to  be  negative  definite  for  all  time. 

\Bl(t)[A(t)p-'  +  P-1AT(t)]BJ.(t)  (2.25) 

This  fact  was  established  by  Fields  [15]  and  is  true  for  the  following  reasons.  Condi- 
tion (2.24)  is  true  if  and  only  if 

{P^B1{t)1x)T(PA(t))(P^B±(t)fi) 

{f-lBMriT{p-iBx(t)J)  ~<0V^0,V*  (2.26) 

Since  the  denominator  of  the  above  fraction  is  positive,  with  some  simplification  the 
above  statement  is  equivalent  to 

^T Bl{t)[A{t)p-'  +  P-Ur(*)]£?.L(^  <  0  V  fi  jt  0,  V  t  (2.27) 

which  is  equivalent  to 

^BlitM^P-1  +  P-Mr(0]5±(i)  <  0  V  t  (2.28) 

If  a  constant  positive  definite  P  exists  satisfying  condition  (2.28),  then  the  given  time 
varying  system  can  be  stabilized.  This  result  is  more  general  than  a  similar  one  given 
by  Shahruz  [12],  and  the  two  are  equivalent  when  P  =  I. 

Before  giving  a  proof  of  Theorem  1,  we  give  an  intuitive  explanation  of  Theorem 
1  for  the  second  order  case  with  a  single  input.  Figure  2.7  shows  the  B  vector  for 
such  a  system  along  with  some  ellipsoid  xTPx  =  1.  If  a  state  vector,  x,  is  drawn 
to  some  point  on  this  ellipse  whose  normal  is  not  perpendicular  to  B,  then  there  is 
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always  a  control,  u,  which  can  place  x  inside  the  ellipsoid.  This  is  true  because  the 
velocity  field  can  always  be  controlled  in  the  B  direction.  This  is  not  true,  however, 
for  points  on  the  ellipsoid  whose  normal  is  parallel  to  span(Bx)  as  shown  in  figure 
2.7.  These  points  on  the  ellipsoid  can  be  found  by  computing 

P^spaniBj,)  (2.29) 

where  B±  is  a  basis  of  column  vectors  orthogonal  to  span(B).  We  now  give  the  proof 
of  Theorem  1 . 

Proof: 
Case  1:  x  €  S(t) 

Substituting  2.20  into  xTPx,  we  have 

xTPx  =  xTPA(t)x  +  xTPB(t)u  (2.30) 

Since  x  €  <5  then 

x  =  P-H(5x(t))  (2.31) 

where  l(Bx(t))  is  any  linear  combination  of  BL{t).  So 

xTPB(t)u    =    liB^fp-^PB^u  (2.32) 

=    \(BL(t)fB(t)u  (2.33) 

=    0  (2.34) 

and  therefore  from  2.30  we  have, 

xTPx  =  x.TPA(t)x  (2.35) 
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Case  2:  x  £  $(t) 

Suppose  we  want  to  set  xTPx  =  c  where  c  is  some  arbitrarily  chosen 
value.  Then  we  want 

c  =  xTPA(*)x  +  xTPJ3(0u  (2.36) 

Since  x  £  5(i)  then  xTPB(t)u  /  0  and  there  exists  at  least  one  u  such 
that  2.36  is  true. 

One  could  solve  (2.36)  for  u  to  use  as  a  control  law,  but  this  would  not  be  practical 
to  implement  because  A  and  B  can  not  easily  be  computed.  Linear  state  feedback  is 
more  desirable.  The  next  section  will  develop  a  procedure  for  computing  the  set  of 
feedback  gains  that  will  implement  a  feedback  control  law  keeping  the  normal  velocity 
component,  xTPx,  within  some  specified  limit. 

2i3 A  Linear  Feedback  Set  to  Control  xTPx 

Recall  from  Section  2.1  that  the  stability  of  a  time  varying  linear  system 

*  =  A(t)x  (2.37) 

could  be  analyzed  by  evaluating  the  bounds  of  xTPA(t)x  on  the  spheroid 

xTx  =  1  (2.38) 

Now  we  want  to  find  the  set  of  all  linear  state  feedback  gains  for  the  control  law, 

u  =  -K(t)x  (2.39) 

for  the  system 

x  =  A(i)x  +  B(t)u  (2.40) 

so  that  condition  (2.19)  holds  for  the  closed  loop  system  where  c  and  c  are  now 
specified.  This  section  will  give  conditions  for  the  nonemptiness  of  such  a  set  followed 


24 


by  a  discussion  of  its  convexity.  Finally  this  section  will  present  an  exhaustive  search 
method  for  finding  the  boundary  of  this  set. 

For  the  remainder  of  this  chapter,  A  and  B  will  be  frozen  at  one  instant  in 
time.  A  discussion  follows  on  how  to  find  a  controller  or  set  of  controllers  which 
satisfies  (2.19)  at  one  given  instant.  If  a  given  system  is  time  varying  then  the  results 
which  follow  must  be  applied  for  all  time  with  P  being  constant  and  positive  definite. 
These  results  apply  to  all  systems  which  can  be  stabilized  with  respect  to  a  Lyapunov 
function  using  a  constant  P.  The  success  of  these  methods  depends  on  the  existence 
of  a  constant  positive  definite  P  so  that 

\Bl{t)[A{t)P-^  +  P-lAT(t)]BL{t)  <  0  V  t  (2.41) 

Since  we  require  P  to  be  zero,  these  results  are  conservative. 

Let  K  and  X  be  the  set  of  all  K  which  satisfy  the  left  and  right  inequalities 
respectively  of  the  following  expression. 

c  <  -xT[PA  +  ATP  -  PBK  -  KTBTP]x  <  c  V  x  :  xTx  =  1  (2.42) 

where  P  is  positive  definite.  The  objective  is  to  find 


K,  :=  JCHJC  (2.43) 

The  following  theorem  gives  conditions  for  the  nonemptiness  of  K  and  X 
Theorem  2  Given  the  'nth  order  m-input  system, 

x  =  Ax  +  Bu  (2.44) 

K  and.  K  are  nonempty  if  and  only  if  there  exists  a  positive  definite  P  so 
that  the  following  conditions  are  true: 

c  <  rmnxTP.4x  =  mm^xT[Pi  +  ATP}x  (2.45) 
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c  >  maxxTP/lx  =  max  \xT[PA  +  ArP]x  (2.46) 


Comments: 

Conditions  (2.45)  and  (2.46)  are  respectively  equivalent  to 


and 


where 


Q  <   \rmn[-{^H)-lBl[AP-1  +  P"1^]  B  L  (  V^)"1]  (2.47) 


>  A^t-fv^fj-^JfAP-1  +  P-Mr]5±(v^)-1]  (2.48) 


H=(p-iBLfp-iB±  (2.49) 

The  matrix  #  is  square,  full  rank,  and  has  dimension  n  -  rn.  The  above  is  true  for 
the  following  reasons.  Equation  (2.46)  can  be  rewritten  as 


C     >  max  —  v* 


=p-™>x=12X  ^A  +  A'/Tx  (2.50) 

=    m^Sx  2^^[j4P_1  +  P"MT^^  (2-51) 

Since  if  is  positive  definite  then  VH  exists,  is  square,  and  has  an  inverse.  By  making 
the  substitution,  p  =  [VS)"lx,  (2.51)  becomes 

~°  ~  ztS  ^(^r'BliAP-1  +  P-1A)Bt(VH)-1z  (2.52) 

which  is  equivalent  to 

c>  Amax[-(V/F)-1^[AJP-1  +  P-'AIB^VH)-1}  (2.53) 

Equation  (2.45)  is  equivalent  to  (2.47)  for  similar  reasons.    We  now  give  a  proof  of 
Theorem  2. 
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Proof: 

We  will  show  the  nonemptiness  of  X.  £  can  be  shown  to  be  nonempty  in 
a  similar  manner. 


Case  1: 


then 


maxxTPAx  <  c 
xes  — 


xTPAx  <  c  V  x  e  S. 


(2.54) 


(2.55) 


Let  K  =  kBTP  where  k  is  a  scalar  value.  Then 

xTPx    =    ^xT[PA  +  ATP  -  PBK  -  KTBTP]x  (2.56) 

=    ^xT[PA  +  ATP}x  -  lk[xTPB][BTPx)         (2.57) 
Since  xTPB  =  [BTPx]T  then 

xTPBBTPx>0\/xgS,   VxTx  =  l  (2.58) 


m 
mce 


If  (2.55)  is  true  and  if  k  is  made  large  enough,  then  the  second  term 
(2.57)  will  dominate  V  x  £  S  and  xTPx  <  c    V  x  :  xTx  =  1.    S 
kB  P  g  JC  then  K.  is  nonempty. 
Case  2: 

maxxTP^x  >  c  (2.59) 

Then 

rp 

x   PAx  >  c  for  some  x  g  S.  (2.60) 

By  Theorem  1,  xPAx  can  not  be  controlled  V  x  g  <5;  so  X  is  empty. 
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Remarks: 

The  nonemptiness  of  the  intersection  K  =  £  n  X  is  not  guaranteed.  If  the  designer 
discovers  that  no  intersection  exists,  then  the  upper  and  lower  velocity  bounds  will 
have  to  be  adjusted. 

Another  useful  property  of  K  is  convexity.  This  property  is  valuable  in  formulating 
iterative  search  techniques  to  be  described  in  the  next  section.  fC  is  convex  if  when 
x  and  y  are  elements  of  K,  then  ax  +  (1  -  a)y  is  also  an  element  of  K  for  0  <  a  <  1 

[16]. 

Theorem  3  Let  X  be  the  set  of  all  K  so  that 

Ka,[-[P{A  -  BK)  +  (A  -  BKfP}}  <  c  (2.61) 

Let  K  be  the  set  of  all  K  so  that 

Kin[-[P{A  -  BK)  +  (A-  BKfP]}  >  c  (2.62) 

Then  /C,  £,  and  K,  =  JC  D  X  are  convex. 

Proof: 

If  £  and  £  are  convex,  then  the  intersection  K,  is  convex.  Convexity  of 
/C  will  be  proven  here.  In  a  similar  manner  the  proof  for  the  convexity  of 
K  can  be  written. 

Let  Kx  and  K2  be  elements  of  X.  We  must  show  that  aKx  +  (1  -  a)K2 
is  contained  in  X  when  0  <  a  <  1.  We  know  that 

IxT[iM  +  ATP  -  PBKX  -  K?BTP}x  <  c  V  xTx  =  1  (2.63) 


and 

-x*[PA  +  ATP  -  PBK2  -  K?BTP]x  <  c  V  xTx  =  1  (2.64) 


1  Jrp  a   ,    ^t. 
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Then 

-axT[PA  +  ATP  -  PBK,  ~  K?BTP]x  <  ac   V  xTx  =  1         (2.65) 
and 

^l-a)xT[PA+ATP-PBK2-I^BTP]x  <  (l-a)c  VxTx  =  l  (2.66) 

where  0  <  a  <  1,  Adding  (2.65)  to  (2.66)  gives 

f xT[/M  +  ATP]x  -  laxT[PBK\  +  Kl BTP]x 

~l(l-a)xT[PBK2  +  iqBTP]x   <c  V  xTx  =  1 
which  becomes 

^xT[PA+ATP-PB(aK1+(l-a)K2)-(aKl+(l-a)K2)TBTP]<c^xTx  =  l 

(2.68) 
So  aKi  +  (1  -  a)K2  €  X  for  0  <  a  <  1.  X  is  convex. 

Now  that  we  have  conditions  on  the  nonemptiness  and  convexity  of  K  and  K  it 
would  be  desirable  to  find  the  boundary  of  these  sets.  We  know  that  if  one  of  the 
eigenvalues  of  the  square  matrix  Q  is  c  then 

det[Q  -  d]  =  0.  (2.69) 

This  is  helpful  in  understanding  the  following  exhaustive  search  method  for  comput- 
ing the  boundary  of  X.  A  similar  procedure  exists  for  finding  dK. 

1)  All  eigenvalues  of 

\[PA  +  ATP  -  PBK  -  KTBTP]  -  cl  (2.70) 

need  to  be  less  than  or  equal  to  zero.  Fix  all  but  one  of  the  entries  of  K. 
Let  the  i,fth  entry  be  the  free  entry  and  be  represented  by  k.  Let  KQiti 
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contain  all  the  fixed  entries  of  K  and  be  defined  as  follows. 


Then 


where 


K0i,j  :- 


ki 


1  K\j 


kit 


kl  0  kin 


"'ml  ""77ij  i^mn 


K  =  kQij  +  Koij 


1        j        n 


Qi,j  '■  — 


m 


0 

0 
0 

0 

0 

1 

0 

0 

0        0        0 


(2.7i; 


(2.72) 


.73) 


2)  Find  all  values  of  the  free  entry  which  make 

det[PA  +  ATP  -  Tel  -  PBK  -  KTBTP]  =  0 


(2.74) 


This  problem  reduces  to  finding  the  roots  of  a  polynomial.  After  sub- 
stituting equation  (2.72)  into  the  argument  of  the  above  determinant,  it 
becomes 

k(-PBQid  -  Ql3BTP)  +  (PA  +  ATP  -  PBKkj  -  X&jEPP)    (2.75) 


Let 


Pi(h3)-=-{PBQUj  +  Ql.BTP) 
Po(Koitj)  :=  PA  +  ATP  -  PBKotj  -  K^B1 'P 


(2.76) 


(2.77) 
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Then  (2.74)  becomes 

detikP,  +  P0]=Q  (2.78) 

Solve  (2.78)  for  k. 

3)  Reject  all  complex  roots.  If  all  the  roots  are  complex  then  skip  the 
next  step. 

4)  Test  the  intervals  between  the  real  roots  by  checking  to  see  if 

Amax[PA  +  ATP  -  Tel  -  PBK  -  KTBTP]  <  0  (2.79) 

The  A"s  that  bound  the  interval  which  satisfies  (2.79)  lie  on  the  boundary 
of  K,.  Convexity  of  /C  implies  that  no  more  than  two  iTs  bound  this 
interval. 

5)  Repeat  the  process  for  all  possible  values  for  the  fixed  entries  in  K. 
The  result  is  dfC. 

dIC  can  be  computed  in  a  similar  way  by  reversing  the  inequalities  in  the  above 
procedure  and  by  replacing  Xmax  with  Xmin.  To  find  K  the  intersection  of  £  and  X 
can  be  found  in  step  4). 

The  fixed  entries  in  step  5)  must  be  assigned  to  a  finite  number  of  grid  points  if 
the  above  procedure  is  to  be  executed  on  a  real  system.  The  spacing  of  these  grid 
points  must  be  smaller  than  the  size  of  the  feedback  set.  If  the  system  has  a  high 
order  or  a  multiple  number  of  inputs,  the  number  of  grid  points  will  become  too  large, 
and  it  will  not  be  practical  to  implement  this  method.  Understanding  this  procedure, 
however,  leads  to  the  formulation  of  an  iterative  method  that  can  be  used  on  high 
order,  multiple  input  systems  and  will  be  presented  in  the  following  section.  First. 
an  example  is  given  applying  the  exhaustive  search  method  to  a  second  order,  single 
input  system. 
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Example  1 

Given  the  system 

where 
and 

let 


x  =  Ax  +  Bu 


(2.80) 


A  = 


B 


-10      12 
0    -1 


cos (75°) 
sin(75°) 


P  =  I 


We  want  to  find  all  feedback  gains  which  satisfy  the  following  constraint. 


c  <  xrPx  <  c  V  x  :  xTx  =  1 


(2.81) 


Before  specifying  c  and  c,  we  must  check  the  normal  velocity  components 
for  x  E  <5  as  Theorem  1  requires.  Let 


x,=Bj  = 


■  sin  75° 
cos  75° 


Then 


xfAxi  =  -12.40 
From  Theorem  2  we  must  have 


(2.82) 


(2.83) 


c<  -12.40  <c 


From  this  we  choose 


(2.84) 


-14 
-9 


(2.85) 
(2.86) 
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In  this  example  n  =  2  and  m  =  1.  It  was  decided  to  set  %  =  j  =  1  so 


that 


Qij  —  Qi,i  — 


1    o 


Koi,j  =  #01,1  =     0    #2 


and 


#  =  kQ1?1  +  K0lil  = 

From  equations  (2.76)  and  (2.77)  come 


k    K2 


Pic 
Pic 

Pqc 
Pqc 


-[BQi,i  +  QixBT] 

-[BQ^  +  Ql.B7] 

A  +  AT  -  2cl  -  BK01tl  -  K*U1BT 

A  +  AT-  Tel  -  BK01A  -  K&  x£r 


(2.87) 
(2.88) 

(2.89) 

(2.90) 

(2.91) 
(2.92) 
(2.93) 


The  roots  of  the  following  polynomials  are  computed  in  terms  of  k  while 
incrementing  K2  through  a  wide  range  of  values. 


det[kPls  +  Poc]    =    0 
detlkP^  +  Pas]    =    0 


(2.94) 
(2.95) 


Rejecting  complex  roots  and  checking  the  regions  separated  by  the  real 
roots  give 


k,(A2),    k£(A2) 
ife(A'2),    %z(K2) 

The  intersection  of  these  regions  are  found. 

k(A2)     =    rnaxlkcikc 
k(A'2)     =     mmfkc,  ]&~ 


(2.96) 
(2.97) 


(2.98) 
(2.99) 
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2       11 


Figure  2.8.  A  plot  of  the  boundary  of  AC  which  guarantees  satisfaction  of  the  design 
constraints  of  Example  1. 

Finally, 


dfC  =  {\k(K2),K2]  VK2}u{[k(K2),K2]  V K2]  (2.100) 

A  plot  of  dfC  is  shown  in  Figure  2.8.  Let  K  =  [7  11].  We  can  check  to  see 
if  K  6  AC  by  evaluating 

Ai    =    \min(A  +  AT  -  BK  -  KTBT)  (2.101) 

=    -12.91  (2.102) 


A2    =    Xmax(A  +  AT-BK~KTBT) 
=    -10.52 


As  the  following  shows, 


(2.103) 

(2.104) 


c    <    Ai 

A2    <    c 


(2.105) 
(2.106) 


34 

2A One  Linear  Feedback  Matrix  to  Control  xTfx 

The  search  procedure  given  in  section  2.3  becomes  impractical  for  high  order 
systems  with  multiple  inputs  because  the  number  of  grid  points  for  the  fixed  entries 
of  K  becomes  very  large.  This  section  discusses  an  iterative  Lyapunov  design  method 
which  saves  on  computations  and  finds  one  element  K  <E  £,  if  it  exists,  where  K  = 
Kl  n  X.  By  applying  this  procedure  at  every  instant  in  time  to  a  time  varying  linear 
system,  one  can  find  a  control  law  stabilizing  the  system  if  a  constant  positive  definite 
P  exists  which  satisfying  the  following  condition. 

\Bl{t)[A{t)P-1  +  P-lAT(t)]BL(t)  <  0  V  t  (2.107) 

By  Theorem  2,  this  condition  guarantees  the  existence  of  a  c  <  0  such  that  X(t) 
is  nonempty  for  every  instant  in  time.  This  procedure  applies  to  all  systems  which 
can  be  stabilized  with  respect  to  a  Lyapunov  function  given  by  a  constant  positive 
definite  P.  The  following  is  a  discussion  of  the  iterative  Lyapunov  method  followed 
by  the  algorithm  itself.  An  example  is  then  given  applying  this  procedure  to  one 
operating  point  of  a  fifth  order  linearized  model  of  the  EMRAAT  missile. 

Figure  2.9  illustrates  the  iterative  procedure.  JCf  is  the  feedback  set  which  satisfies 
the  designer's  predetermined  constraints  for  some  specified  P.  The  constraints  are 

cf  <  xT[PA  +  ATP  -  PBK  -  KTBTP]x  <  cf  V  x  :  xTx  =  1  (2.108) 

Let  fd  be  defined  as  the  set  of  all  K  which  satisfies  the  following  constraints. 

Si  <  xT[PA  +  ATP  -  PBK  -  KTBTP]x  <  c{  V  x  :  xTx  -  1  (2.109) 

Given  Kh  we  would  like  to  find  q  and  c{  so  that  if  K%  $  £,  then  K%  e  did.  We  also 
require  that  Kf  C  Kx.  If  K%  £  JCf,  then  we  want  £,■  =  £/.  The  following  definitions 
for  Ci  and  ct  meet  these  requirements.  Let 

*•'  :=  iSS  l*Tip(A  -  BKi)  +  {A-  BKifPjx  (2.110) 

11*1 1  —  1  L 
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Figure  2.9.  A  geometric  view  of  the  iterative  Lyapunov  design  method. 
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and 

k  :=  iSi=i  \*TtP(A  -  BKt)  +  (A  -  M,)2* P]x  (2.111) 

Then  let 

Ci  =  max[cf,\i]  (2.112) 

and 

fii  =  min\cf,  A,-]  (2.113) 

A'0  is  the  initial  guess  in  the  search  for  Kf  6  Kf.  c0  and  Cq  are  computed  using 
(2.112)  and  (2.113)  so  that  K0  is  a  member  of  8KQ  and  £/  C  KQ.  Then  a  new 
feedback  matrix,  Kx  is  found  which  lies  inside  of  JC0,  but  not  on  the  boundary.  New 
constraining  values  are  found  in  the  same  way  as  before  so  that  the  boundary  of  the 
next  feedback  set  contains  Kx.  K2  is  then  found  so  that  it  lies  inside  of  the  present 
feedback  set,  but  not  on  its  boundary.  This  process  is  continued  until  K{  =  Kf  £  fCf. 
The  success  of  finding  Kf  depends  on  the  following  conditions. 

1.  Given  that  K{  e  did  we  must  be  able  to  find  Kl+1  such  that  K,+l  £  AC,-. 

2.  We  must  show  that  ci+1  <  c,  when  c%  >  c/  and  ci+1  >  c{  when  c{  <  cf. 

3.  We  must  show  that  Kf  C  JCi+1  C  K.{. 

4.  Kf  must  be  nonempty. 

We  now  address  these  four  points. 

1.  Given  Ki  e  dKi,  we  need  to  find  a  second  feedback  matrix  Kl2  £  8Kt  where 
Ki2  +  Ki,  Then,  due  to  convexity,  \K%  +  \K%2  is  a  member  of  Kt.  Figure  2.10a  shows 
a  second  order  example  of  a  procedure  for  finding  Kl+1.  The  algorithm  will  be  given 
shortly.  The  horizontal  and  vertical  axis  are  assigned  to  kn  and  kl2  respectively.  The 
region  Kt  is  enclosed  by  dK,  and  d~Kx.  A't  is  known.  Kt2  is  found  by  searching  along 
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the  line  that  passes  through  Kt  and  is  parallel  to  the  kn  axis.  Point  c  is  found  by 
computing  \Ki  +  \K%2.  Since  lCt  is  convex,  then  c  €  K,%.  This  step  is  repeated  again 
by  searching  along  the  k12  axis.  Using  similar  arguments,  point  /  is  also  in  K%.  Ki+1 
is  set  equal  to  point  /.  For  higher  order  systems,  additional  boundary  points  are 
found  and  interior  points  computed  by  searching  along  directions  which  are  parallel 
to  the  axis  of  the  coordinate  system. 

Under  normal  conditions,  this  procedure  works  well.  However  numerical  problems 
do  occur.  These  problems  will  now  be  discussed  along  with  a  cure.  Figure  2.10b  shows 
a  second  order  example  of  when  the  above  method  fails.  The  boundaries  are  shaped 
m  such  a  way  that  searching  along  line  /,  and  l2  yields  no  new  boundary  points.  A 
solution  to  this  problem  is  to  relax  one  of  the  constraining  values  by  setting,  for  ex- 
ample, c{  equal  to  c{  +  6.  dZ{  will  then  move  so  that  point  b  can  be  found  as  shown  in 
figure  2.10c.  Later,  the  relaxed  constraining  value  can  take  on  its  original  assignment. 

2.  We  need  to  show  that  c,-  >  ci+1  when  c,  >cf.  In  this  case 

Vi=li>Zf  (2.114) 

Since  Ki+1   <S  Id  and  Ki+l   £  did,  then  the  following  condition  holds  with  strict 
inequality. 

-xT[P(A  -  BK,+1)  +  (A-  BKi+1fP]x  <c,Vx:xTx=l  (2.115) 

Taking  (2.110)  for  i  +  1  then 

.  ■     -xT[P(A-BK,+1)  +  (A-BKi+1)TP}x<Jt+1yX:xTx  =  l  (2.116) 

and  there  exists  an  x  which  satisfies  the  above  condition  for  equality.    Therefore 
ct  >  Xi+1  and  C{  >  max[cf,Xi+1]  =  cl+l.    Using  similar  arguments,  it  can  be  shown 
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"12 


dJC,   ^ 


a»- 


IKi 


/  dXi 


I 


:*— - 


Hi 


fell 


Figure  2.10.   Step  5  of  the  iterative  Lyapunov  design  method  for  (a)  a  second  order 
example,  (b)  how  it  sometimes  fails,  and  (c)  how  this  problem  is  corrected. 
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that  c?:  <  cl+l  when  c{  <  C/. 

3.  We  now  show  that  fCf  C  K,l+1  C  £,-.  Since  c,-  =  max[c/,  A,-],  then  ct  >  cf.  We  have 
already  seen  that  q  >  A4+1.  So  ct  >  rnaxfc/,  Xi+l]  =  c!+1.  Similarly  c,-  <  ct+1.  /Ct+1  is 
the  set  of  all  K  so  that 

ci+1  <  ^T[P(A  -  BK)  +  (A  -  BK)TP]x  <  cl+1  V  x  :  xTx  =  1  (2.117) 

Since  Ci  >  ci+1  and  ct  <  ci+1,  then  for  every  element  of  K,l+X  the  following  holds. 

Si  <  \*T[P(A  -  BK)  +  (A-  BKfP]x  <c,Vx:xTx=l  (2.118) 

Therefore,  fCi+1  C  JC{.  Since  ci+1  >  cj  and  ci+1  <  cf,  then  using  a  similar  argument 
fCf  C  fCi+i- 

4.  In  using  this  design  procedure,  P  is  chosen  so  that  the  maximum  uncontrol- 
lable normal  velocity  component  is  negative.  Then  from  Theorem  2,  cs  can  be  made 
negative  in  an  attempt  to  achieve  stability.   cf  must  be  greater  than  the  maximum 
uncontrollable  normal  velocity  component,  and  cf  must  be  less  than  the  minimum 
uncontrollable  normal  velocity  component.  From  Theorem  2  this  will  guarantee  the 
nonemptiness  of  JCj  and  Kf.  The  nonemptiness  of  the  intersection  of  these  two  sets, 
however,  is  unknown.   If  /C/  is  nonempty,  then,  as  i  becomes  large,  K,  €  JC/.  If  JC/ 
is  empty  then  c,  and  c{  will  converge  to  values  which  do  not  match  the  desired  con- 
straints and  Ki  will  yield  a  closed  loop  system  that  meets  the  constraints  given  by  cz 
and  Cj.  The  designer  will  either  have  to  accept  this  result  or  try  again  with  a  different 
P  or  different  constraining  values  or  both.    Since  stability  is  desired,  one  approach 
would  be  to  keep  P  and  cf  and  lower  cf  until  ICf  becomes  large  enough  to  intersect  Kj. 

The  outline  of  the  iterative  Lyapunov  design  method  is  as  follows. 
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1)  Choose  P.  Cf,  and  cf.  This  selection  must  obey  Theorem  2.  It  should 
be  noted  that  Theorem  2  guarantees  the  nonemptiness  of  JCf  and  £/,  but 
not  their  intersections.  If  the  system  is  time  varying,  then,  in  order  to 
use  this  algorithm,  P  must  be  found  so  that 

l-Bl{t)[A{t)p-1  +  p-1AT(t))BL(t)  <  0  V  t  (2.119) 

Otherwise  this  algorithm  cannot  guarantee  stability. 

2)  Compute  c0  and  c0  so  that  A'0  =  0  e  <9/C0  and  K}  C  K,0.  K0  will  be 
the  initial  guess  in  the  search  for  Kf  €  ICj. 

3)  Let  t  =  0 

4)  Let  i  =  i  +  1 

5)  Find  j$  so  that  A",  g  JCt-_!  and  A't  ^  a^_a. 

6)  Compute  ct  and  c,  so  that  A",  €  dK{  and  X!/  C  C»  or  so  that  K%  e  JCf. 

7)  Repeat  steps  4)  5)  and  6)  until  one  of  two  events  occur. 
1.  ct  =  Cf  and  c,-  =  Cf. 

2-  (Cj  —  £s-_x)  and  (e;  —  Cj_i)  become  very  small. 

Remarks: 

If  event  1  occurs,  then  fCf  is  nonempty  and  A,  e  fCf.  If  event  2  occurs,  then  JCy  is 
empty  and  Kz  yields  a  closed  loop  system  that  meets  the  constraints  corresponding 
to  c,-  and  c,-. 

We  now  explain  how  to  perform  steps  2),  5),  and  6). 


Step  2) 
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The  lower  and  upper  bound  of  xTPx  for  the  open  loop  system  is  resp 
tively, 


ec- 


Then, 


Ao  =  Amm[-(PA  +  ArP)] 
A0  =  AmQX[-(PA  +  ArP)] 


(2.120) 
(2.121) 


c0  =  mmfc^Ao] 
c0  =  max[~Cf,  X0] 


(2.122) 
(2.123) 


Step  5) 


The  algorithm  is  now  given. 


K  =  Ki 

For  j  =  1  to  m 

For  k  =  1  to  n 

K0j,k  =  K  -  kJikQjtk 

Pi  =  -{PBQhk  +  QlkBTP) 

Po  =  P(A  -  BK0j,k)  +  {A-  BKQhk)TP  -  2cl 

Solve  det[kPx  +  P0]  =  0  in  terms  of  k  for  c  =  ci?  and  c{. 

Reject  all  values  which  do  not  meet  the  constraints  from 

steps  2)  and  6).  The  result  is  two  intervals  whose 

lower  bound  is  k±  and  h^  and  whose  upper  bound  is 

k\  and  hi- 

Find  the  intersection  of  these  intervals  by  evaluating 
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k  =  max[kuk^[ 
k  =  min[ki,  k2]. 
Find  the  midpoint  of  the  interval  by  computing 

Let  A'  =  kj,kQj,k  +  Kojtk 
Next  fe 
Next  j 
Ki+i  =  K 


Step  6) 


Let 


A    =    \rmn[\(PA  +  ATP-PBKt-K?BTP)}  (2.124) 

A    =    A^f-CPA  +  ^P-Pfiifi-iirf^P)]  (2.125) 


So, 


Si  =  mm[£/5A]  (2.126) 

Ci  —  max[cf,X]  (2.127) 

The  following  example  illustrates  the  feasibility  in  applying  this  method  to  the 
EMRAAT  missile. 

Example  2 

We  would  like  to  apply  the  iterative  Lyapunov  design  method  to  the  EM- 
RAAT missile.  The  missile  was  flown  in  a  simulation  through  a  trajectory 
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using  another  autopilot  design.    The  model  was  linearized,  and  the  fol- 
lowing system  was  taken  at  4.00  seconds  into  the  flight. 


x  =  Ax  +  Bu 


(2.128) 


where 


A  = 


-.9939 
.1613 
-64.36 
-252.2 
-.5806 


-.1612 

-.4100 

-1557 

-6.097 

83.29 


.0086 
.1079 

-2.078 
.0173 

-.0261 


.9997 
0.0 
.1763 
.6403 
.1597 


0.0 

-.9997 

.0442 

.1611 

-.5701 


(2.129) 


B  = 


0.0 

-.0149 

-1150 

-4.494 

1.234 


.1296 

0.0 

0.0 

.1044 

31.06 

-1222 

121.7 

-4.747 

.2802 

-108.6 

(2.130) 


and 


x  =  [a,f3,p,q,r}' 


(2.i3r 


Step  1) 
Let  P  =  I.  The  minimum  and  maximum  uncontrollable  normal  velocity 
components  were  found  by  computing  the  limiting  values  in  equations 
(2.47)  and  (2.48).  Since  P  =  /,  these  terms  simplify  and  are  evaluated  as 
follows: 


^m-xT[PA  +  ATP]x 


^rmn[-Bl(A  +  AT)B±]       (2.132) 
-•7263  (2.133) 


and 


m&x-xT[PA  +  ATP]x 
xes  2 


KaX[-Bl(A  +  AT)BL]       (2.134) 
-•3017  (2.135) 
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Theorem  2  guarantees  the  nonemptyness  of  K.  and  fC  if  cf  >  -.3017  and 
cj  <  -.7263.  With  this  in  mind,  we  let  cj  =  -.3  and  cf  =  -6000.  Since 
the  intersection  of  K.  and  K  is  not  guaranteed,  cf  was  chosen  to  be  very 
negative  to  increase  the  probability  of  getting  an  answer. 
Step  2) 


A      =      Xrmn[-(A  +  AT)]  (2.136) 


•2 
=    -781.4  (2.137) 


So 


"1 


138) 


A     =     \max[^(A  +  AT)}  (2.139) 

=    778.9  (2.140) 

(2.141) 


c0     =    min(cf,X)  (2.142) 

=    -6000  (2.143) 

c0    =    max(cf,X)  (2.144) 

=    778.9  (2.145) 

Steps  3)  -  7) 

A  program  was  written  for  MATLAB  to  carry  out  the  iterations  in  steps 
3)  -  7).  The  program  would  terminate  if  fCj  was  found  or  when 

8  :=  \Zi  -Zi-x\  <  t  =  .0001  (2.146) 


JCf  was  found  before  6  <  .0001.  The  program  ran  16  iterations.  The  final 
result  is 


Kf  = 


.00046       2.125    -.7974       .0805       .7535 

2.057    -.0012       2.339    -6.577       .2591 

-.00029    -.7259       .0263    -.1346    -.7985 


Now  to  check  the  result. 


(2.147) 
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A,      =      \mink(A  +  AT-BKf-KJBT)} 


=    -818.0 


(2.148) 
(2.149) 


A/    =    Xmax[-(A  +  AT-BKf-KJBT)} 


=    -.3001 


(2.150) 
(2.151) 


This  meets  the  desired  constraint. 


Qf  <  Xf  <  Xf  <  cf 


CHAPTFR  3 
A  TIME  VARYING  SECOND  ORDER  EXAMPLE 

The  following  problem  gives  a  case  when  pole  placement  succeeds  in  giving  eigen- 
values with  negative  real  parts,  but  fails  to  stabilize  the  system.  The  Lyapunov 
design  method  is  then  employed,  and  the  resulting  closed  loop  system  is  shown  to  be 
asymptotically  stable  for  all  time. 

We  would  like  to  find  a  feedback  control  law  that  stabilizes  the  systt 


tern 


x  =  A(t)x  +  B(t)u 


(3.1) 


where 


A(t)  = 


and 


-1  +  1.5cos(<)[cos(*)  +  cos(<  +  tt/18)]       1  -  1.5sin(t)[cos(i)  +  cos(t  +  tt/18)] 
-1  -  1.5cos(i)[sin(i)  +  sin(f  +  tt/18)]    -1  +  1.5sin(t)[sin(i)  +  sm(t  +  tt/18)] 

(3.2) 


B(t)  = 


cos{t  +  tt/18)  1 
■  $m(t  +  tt/18)  I 


(.3.3) 


The  eigenvalues  of  A(t)  are  -.4889  and  1.4661  for  all  time.  The  following  control  law 
is  proposed. 

u  =  [  L5cos(i)    -1.5sin(t)  ]x  (3.4) 

The  resulting  closed  loop  system  can  be  found  in  example  5.3,109  by  Vidyasagar 
[10]  and  also  in  Khalil  [13].  The  eigenvalues  for  the  resulting  closed  loop  system  are 
A  =  -.25  ±  j.6614.  Since  the  eigenvalues  have  negative  real  parts,  one  would  expect 
the  closed  loop  system  to  be  stable.  However,  Vidyasagar  shows  that  the  transition 
matrix  is 


$(t,0)  = 


e-5icos(i)     e-*sin(^) 
-e-5tsin(i)     e-*cos(t) 


(3.5) 
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300 


Trajectory  of  a  pole  placement  time  varying  system 


600 


Figure  3.1.  The  trajectory  of  the  above  closed  loop  system  based  on  pole  placement. 
The  initial  conditions  are  x0  =  [1  0]r.  The  system  is  unstable. 


is 


If  the  initial  conditions  of  the  system  are  x0  =  [1  0]r,  the  resulting  trajectory  i 
unstable  as  figure  3.1  shows. 

We  now  turn  to  the  Lyapunov  design  method.  We  choose  P  to  be  the  identity 
matrix.  Before  giving  the  design  constraints,  we  need  to  check  the  the  value  of  the 
uncontrollable  normal  velocity  components  for  all  time.  At  t  =  0 


B  =  [.9848  -  .1736]' 


(3.6) 


Since  P  is  the  identity  matrix,  we  are  interested  in  the  normal  velocity  component 
which  is  on  the  part  of  the  unit  circle  whose  tangent  is  parallel  to  B.  So  we  let 


x  =  B±  =  [.1736  .9848]T 
The  uncontrollable  normal  velocity  component  at  t  —  0  is 

xTPx  =  -xr[PA(0)  +  AT(0)P]x  =  -.9548 


(3.7; 


[3.8) 
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The  uncontrollable  velocity  for  this  example  is  constant  for  all  time.  In  selecting  the 
constraining  values,  we  must  have 

c<-.9548<c  (3.9) 

The  following  assignments  are  made 


c=  -2 


(3.10) 


c  =  -.5  (3.11) 

Since  the  system  is  second  order  and  has  only  one  input  it  is  possible  to  plot  the  set 
of  all  feedback  gains  which  satisfy  the  following. 

c<lxT[P(A(t)-B(t)K(t))  +  (A(t)-B(t)K(t))TP}x<c\/x:xTx=l:Vt  (3.12) 

This  time  varying  feedback  set  is  shown  in  figure  3.2.  Since  the  time  varying  nature 
of  the  system  is  periodic,  and  from  inspection  of  figure  3.2,  the  following  control  law 
is  chosen. 

u  =  3[cos(*)   -sin(t)Jx  (3.13) 

The  feedback  matrix  in  (3.13)  is  shown  to  be  inside  the  moving  feedback  set  in  figure 
3.2.  The  eigenvalues  of  the  derivative  of  the  resulting  Lyapunov  function  is 

X[A(t)  -  B(t)K(t)  +  AT(t)  -  KT(t)BT(t)]  =  -1.1193,  -0.8579  Vi  (3.14) 

Therefore,  the  closed  loop  system  is  stable.  Figure  3.3  shows  a  trajectory  of  the 
Lyapunov  based  closed  loop  system  where  the  initial  condition  is  x0  =  [1  0]T. 
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Figure  3.2.  The  time  varying  feedback  set  which  satisfies  the  design 


gn  constraints. 
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0.05 


Trajectory  of  the  Lyapunov  based  closed  loop  time  varying  system 


Figure  3.3.  A  trajectory  of  the  closed  loop  system  using  the  Lyapunov  design  method. 
x0  =  [1  0]r.  The  system  is  stable. 


CHAPTER  4 
DERIVATION  OF  A  MODEL  OF  THE  EMRAAT  MISSILE 

The  next  section  derives  a  nonlinear  model  of  the  EMRAAT  missile.  A  linearized 
model  is  then  generated  in  the  following  section  and  will  be  used  for  the  design  of  the 
autopdot.  Aerodynamic  and  inertial  cross  coupling  are  assumed  negligible  in  order 
to  reduce  the  order  of  the  model.  A  linearized  pitch  model  and  a  linearized  roll-yaw 
model  results.  All  assumptions  are  clearly  stated. 

4.1     The  Nonlinear  Model 

Before  deriving  the  nonlinear  model  some  variables  and  terms  are  defined.  Figure 
4.1  shows  the  missile  body  coordinate  frame  of  the  EMRAAT  missile  [2],  The  three 
axes,  x,  y,  and  z,  are  fixed  to  the  missile  as  shown.  The  velocity  of  the  missile  is 
represented  by  V.  Angle  of  attack,  a,  is  defined  as  the  angle  between  x  and  the 
projection  of  V  onto  the  x-z  plane.  Sideslip,  /?,  is  the  angle  between  x  and  the 
projection  of  V  onto  the  x-y  plane.  The  velocity  components,  u,  v,  and  w,  are  the 
projections  of  V  onto  the  x,  y,  and  z  axis  respectively.  The  angle  rates,  p,  q,  and  r, 
are  named  roll  rate,  pitch  rate,  and  yaw  rate  respectively  and  are  defined  as  the  rate 
of  rotation  around  the  x,  y,  and  z  axis.  Their  directions  obey  the  right  hand  rule. 
These  definitions  are  similar  to  those  applied  to  aircraft  as  given  by  Etkin  [17]. 

The  following  derivation  of  the  nonlinear  model  is  based  on  a  similar  model  found 
m  Smith  [2]  and  in  Koenig  [18].  Assumptions  are  made  here  to  separate  the  system 
into  two  lower  order  models:  the  pitch  model  and  the  roll-yaw  model.  We  begin  the 
derivation  by  starting  with  the  force  equations.  The  forces  acting  on  the  missile  are 
thrust,  gravity,  and  aerodynamic  forces.  The  autopilot  design  in  this  paper  applies 
to  the  second  phase  of  the  flight  when  the  engine  is  no  longer  burning  so  thrust  is 

51 


52 


Figure  4.1.  The  missile  body  coordinate  frame  of  the  EMRAAT  missile  [2]. 
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zero.  Also,  the  weight  of  the  missile  is  small  in  comparison  to  the  aerodynamic  force 
so  gravity  is  neglected.  The  aerodynamic  force  in  the  x  direction  is  much  smaller 
than  the  aerodynamic  force  in  the  y  and  z  directions,  and  therefore,  will  be  ignored 
for  the  rest  of  the  paper.  Newton's  second  law  of  motion  implies  the  following: 


Fy  =  m[v  +  ru  —  pw] 


•1) 


Fz  =  m[w  +  pv  —  qu]  (4.2) 

FY  and  Fz  are  the  aerodynamic  forces  in  the  y  and  z  directions.  The  quantities 
inside  the  brackets  are  the  total  accelerations  in  the  y  and  z  directions.  We  know 
that 

V  =  [u2  +  v2  +  w2]2  (4.3) 

Since  v  and  w  are  much  smaller  than  u  then 


V  ~ 


u 


Angle  of  attack  and  sideslip  are  given  by 


,  w 
a  =  arctan  (  — 


(3  =  arctan  (  — 
u 


If  we  assume  that  a  and  /?  are  small,  then 


and 


w 

a  ~  — 

u 


(3cv- 
u 


Equations  (4.1)  and  (4.2)  can  be  rewritten  as 


Fy  =  mu 


V 

w' 

-+T 

-  p— 

u 

u . 

(4.4) 


(4.5) 

(4.6) 


(4.7) 


(4.81 


(4.9] 
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r>  W  V 

rz  =  rnu f-  p q 

.u         u 

which  simplifies  to 

/->       I,  r^ 

rY  =  mv    -  +  r  —  pa 

■  U 

\w 

Fz  =  mV   —  +p/3  -q 

.u 

We  assume  that  the  forward  velocity  changes  slowly  so  that 


u  ~  0. 


Then 


and 


w 


a  ~ 


u 


u 


So  (4.11)  and  (4.12)  become 


mV 
Fz_ 

mV 


=  $  +  r  —  pa 
—  a  +  p/3  —  q 


The  aerodynamic  forces  are  given  by 


Fy  =  QS  [Cy0/3  +  CYpP  +  CYrr  +  CYJP  +  CYJr 

Fz  =  -QS  [CNaa  +  CNaa  +  CNqq  +  CNfq6q 
Q  is  the  dynamic  pressure  and  is  defined  as 

Q  =  \pV* 


(4.10) 

(4.11) 
(4.12) 


(4.13) 

(4.14) 

(4.15) 

(4.16) 
(4.17) 


(4.18) 


(4.19) 


(4.201 


where  p  is  the  air  density.  S  is  the  surface  area  of  the  wing.  The  aerodynamic 
coefficients  come  from  wind  tunnel  tests.  They  depend  on  mach  number,  and  some 
depend  on  angle  of  attack.  The  values  of  these  coefficients  have  been  put  in  tabular 
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form  and  are  given  in  the  first  appendix  along  with  other  data  pertaining  to  the 
EMRAAT  missile. 

Substituting  (4.18)  and  (4.19)  into  (4.16)  and  (4.17)  and  solving  for  a  and  $  gives 


and 


A   ,    QSCN  QS 

•  ..—a  =  q  —  pp 


a  = 


1  + 


mV 
QSCn6 


mV 


-l 


mV 
QS 


CNaa  +  CNqq  +  CNsSt 


6qUq 


q~P(*~mV  \CN«a  +  CN*q  +  Cn^ 


<5<J      9 


(4.21) 


(4.22; 


(4.23) 


Equations  (4.21)  and  (4.23)  are  two  nonlinear  state  equations.  In  order  to  separate 
pitch  dynamics  from  roll-yaw  dynamics,  0  is  assumed  close  to  zero  in  (4.23).  Thus, 

QS 


a  = 


1  + 


QSCn~  "" 


rnV 


QS;CNa)  a+(l-^CNq)q+  r—,NsqJ  Uq 


mV 


:Cf 


(4.24) 


Equation  (4.21)  can  be  written  as 

*  -  (&*0  »(•  +  $*) »(-'  +  &*)  <3^)  *(&*)  * 

(4.25) 

We  now  turn  to  the  moment  equations  of  the  missile.  Since  thrust  is  zero  and  since 

gravity  does  not  contribute  any  moment  to  the  missile,  then  the  moments  around 

the  x,  y,  and  z  axis  are  given  by  /,  m,  and  n  respectively,  the  moments  due  to 

aerodynamic  pressure.  They  are  given  by 

(4.26) 

(4.27) 
and 


/  =  QSd  \Ch0  +  ClpP  +  Clrr  +  ClsJp  +  CuA 


Sf ,P      <      y"'tSr'JT 


m  =  QSd  [Cmaa  +  Cm,a  +  Cma  +  CmJt 


Sq      1 


n  =  QSd  [Cns/3  +  Cnpp  +  Cnrr  +  CnJp  +  CnsSr 


(4.28) 
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where  d  is  the  missile  diameter.  Again,  the  aerodynamic  coefficients  come  from 
wind  tunnel  testing  and  are  tabulated  in  the  first  appendix.  Euler's  three  moment 
equations  can  be  written  as  follows. 


/  ' 

'  P  ' 

m 

=  J 

q 

n 

r 

+  H(p,q,r) 


(4.29) 


where 


J  = 


Ixx      —IxY 
—  IxY  IyY 

—Ixz    —Iyz 


-Ixz 

-Iyz 

Izz 


(4.30) 


and 


H  = 


(4.3i; 


-(Iyy  -  hz)qr  +  IyZ(r2  -  q2)  -  Ixvpq  +  Ixyrp 
-{Izz  ~  Ixx)rp  +  Ixz(p2  ~  r2)  -  IXyqr  +  IYZpq 
-(Ixx  -  hY)pq  +  Ixv{q2  -  p2)  -  hzrp  +  Ixzqr  _ 
We  will  assume  that  the  inertial  cross  products,  IXy,  Ixz,  and  IYZ  are  small.  Also 

due  to  symmetry  of  the  airframe,  Iyy  and  Izz  are  assumed  to  be  equal.  Inertial  data 

for  the  EMRAAT  missile  is  given  in  the  second  appendix.    Solving  (4.29)  for  p,  q, 

and  r  gives 


p  = 


1 


xx 


-I 


q  = 


i- 


YY 


[m  +  (Izz  -IXx)rp] 


and 


1 


r  — 


I 


zz 


[n  +  (Ixx  -  Iyy)pq] 


(4.32) 
(4.33) 

(4.34) 


Substituting  (4.26),  (4.27),  and  (4.28)  into  (4.32),  (4.33),  and  (4.34)  gives 

QSd 


and 


r  — 


q  = 


QSd 


p  = 


QSd 


I 


xx 


CipP  +  ClpP  +  Clrr  +  ClspSp  +  Clgr6r 


Iyy 


Cmaa  +  Cmaa  +  Cmqq  +  Cms  Sq 


.  Izz  -  Ixx 

— r —  — rp 


I 


YY 


(4.35) 
(4.36) 


^j        [Cnfl/3  +  Cnpp  +  Cnrr  +  Cnsp6p  +  CnSrST    +  -  ^—  —  pq 


I 


zz 


(4.31 
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If  gyroscope  effects  are  considered  small,  then  (4.24),  (4.25),  (4.35),  (4.36),  and  (4.37) 
can  be  written  as  two  separate  systems:  the  pitch  dynamics,  and  the  roll-yaw  dy- 
namics. The  pitch  dynamics  are  as  follows. 


a  — 


1  + 


QSC 


/Vn 


mV 


-l 


QS 
mV 


C 


Na 


a+     1 


QS 
mV 


CN.)q+(-^CNiq)69 


1  = 


QSd 


lYY 


Cmaa  +  Cmaa  +  C    q  +  Cm,S 


6q^1 


(4.38) 
(4.39) 


The  roll-yaw  dynamic  equations  are 


»  =  (£*.)  M-  +  gCr)  »(-l  +  3&A  r+mCYlr)  SA^C,  ]  S, 


QSd 


mV 


\mV 


P 


r  = 


Ixx 
QSd 
I 


(4.40) 
(4.41) 


zz 


?l,0  +  CipP  +  Clrr  +  CltpSp  +  QJr 
Cn0{3  +  CnpP  +  Cnrr  +  CnspSp  +  CntrSr]  (4.42) 

The  states  of  the  pitch  model  are  a  and  q  with  Sq  as  its  input.    For  the  roll-yaw 
model,  the  states  are  /?,  p,  and  r,  and  the  inputs  are  Sp  and  Sr. 

4.2     The  Linear  Model 

The  previous  section  gave  a  model  of  the  pitch  dynamics  and  the  roll- yaw  dynam- 
ics in  the  following  form. 

[  a     q  }T  =  fq(a,q,6q)  (4.43) 

[0    p    r  }T  =  fpr([3,p,r,6p,6r)  (4.44) 

We  would  like  to  have  two  linearized  models  for  use  with  the  proposed  design  method. 
The  result  will  be  two  systems  in  the  following  form. 


x  =  Ax  +  Bu 


(4.45) 
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where  x  contains  the  states  and  u  contains  the  inputs.  A  and  B  are  matrices  which 
are  functions  of  several  time  varying  flight  parameters  and  are  computed  as  follows. 


dx 


(x,w)r 


B 


du 


(4.46) 


(4.47) 


\*i  ™  )  nominaL 

where  w  contains  additional  flight  parameters.   Note  that  x  and  u  are  now  pertur- 
bations from  the  point  around  which  the  linearization  is  taken. 

We  linearize  the  pitch  model  first.  From  inspection  of  (4.38)  we  see  that 

V  mv     J       \  mV 


aqi2  =  (1  + 
bqii  =  ( 1  + 


QSCn6 
mV 

QSCNa^-1 


QSr 
-QSr 


mV    J       \mV     "Sq, 
To  linearize  (4.39),  we  must  first  substitute  (4.38)  in  for  a.  Then  we  differentiate 
in  (4.46)  and  (4.47). 


as 


aq2\  — 


aq22  — 


bq2l  = 


QSd 

lyy 

QSd  ' 

Iyy 

QSd 
Iyy 


^m„  +  ^ma  (  1 
Cma  +  Cm^    (  1  + 


QSCnj, 
mV 

QSCNa 
mV 


-QS 
mV 


C 


Nc 


i     QSr 

1  ~  mV°N< 


Cmf.„  +  Cm^    1  -\ -~        (  — — CV 

mV    J      \  mV       6q 


*6q 


The  resulting  linearized  pitch  model  is. 


a 


aqll      0-ql2 
aq21       dq22 


q 


+ 


bqll 
bq21 


(4.48) 


The  same  procedure  is  applied  to  the  roll-yaw  model.  Assuming  that  a  is  constant 
in  (4.40),  and  from  inspection  of  equations  (4.40)  to  (4.42)  the  following  results. 


^prll   — 


QSCYs 

r      i    <lprl2  —  OL  + 


rn 


QSCy 

mV 
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CHAPTER  5 
THE  DEPENDENCE  OF  GAINS  ON  FLIGHT  PARAMETERS 

This  chapter  describes  the  procedure  that  was  used  to  determine  which  flight 
parameters  would  be  scheduled  against  the  gains.   Figure  5.1  gives  a  block  diagram 
of  the  system  used  for  this  procedure.    The  feedback  gains  are  computed  using  the 
iterative  Lyapunov  design  method  described  in  Chapter  2.    The  gains  depend  on 
the  linear  model  which  in  turn  depends  on  seven  flight  parameters.    Six  of  these 
parameters  are  held  constant  while  the  seventh  one  changes.  The  resulting  gains  are 
checked  to  see  if  they  depend  on  this  changing  parameter.  The  process  is  repeated 
for  each  flight  parameter.  It  is  found  that  the  gains  depend  on  angle  of  attack,  mach 
number,  and  dynamic  pressure.  This  information  will  be  used  in  the  gain  scheduling 
process.  The  next  section  will  discuss  the  atmospheric  tables  used  to  generate  p  and 
V  from  M  and  Q.  The  flight  parameter  generator  will  then  be  presented.  The  third 
section  describes  the  initializer.  The  details  of  the  iterative  Lyapunov  design  method 
are  given  in  the  fourth  section.    Finally,  the  results  of  the  comparison  between  the 
gains  and  flight  parameters  will  be  given  in  the  fifth  section. 

5T Generating  o  and  V  from  M  and  Q 

By  inspection,  the  linear  models  from  Chapter  4  clearly  depend  on  a,  0,  p,  q,  r. 
p,  V,  and  the  aerodynamic  coefficients.  The  coefficients,  however,  can  be  eliminated 
from  the  list  because  they  depend  on  mach  number  and  angle  of  attack.  It  is  desirable 
to  replace  p  and  V  with  M  and  Q  since  the  later  two  can  be  easily  measured  on  the 
missile.  We  know  the  following. 

Q  =  \pV2  (5.1) 

*,   ' ' v 

M  =  —  (5-2) 

v  SOS 
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Figure  5.1.    A  block  diagram  of  the  systemy  used  to  determine  the  dependence  of 
gains  on  flight  parameters. 
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where  Vsos  is  the  speed  of  sound.  Both  p  and  Vsos  are  functions  of  altitude. 

P  =  fP(h)  (5.3) 

Vsoa  =  f,(h)  (5.4) 

where  h  is  altitude  in  feet  above  sea  level.  Here,  fp  and  /,  are  functions  based  on 
atmospheric  tables  and  are  implemented  by  linear  interpolation.  Solving  (5.2)  for  V 
and  substituting  the  result  into  (5.1)  gives 

Q  =  \pM2VL  (5.5) 

We  generate  a  third  table  in  the  following  way. 

Mh)  :=  fP(h)ps(h)  =  PVls  (5.6) 

So 

Q  =  l-M2f,{h)  (5.7) 

The  function,  /3,  is  a  one-to-one  function  so  that  its  inverse  can  be  found  by  reading 
the  table  backwards.  With  this  in  mind  we  can  solve  (5.7)  for  h. 

h  =  k1  iw)  ^ 

From  (5.2)  and  (5.4) 

V  =  MVS0S  =  Mf,(h)  (5.9) 

Substituting  (5.8)  into  (5.9)  to  eliminate  h  gives 

v  m  Mf-  (^  (§))  («-10) 


Also  from  (5.3) 


/>  =  /,( /a-1  (§))  fo.n: 
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Equations  (5.10)  and  (5.11)  are  used  to  eliminate  p  and  V  in  the  linear  model.  As 
a  result,  the  linearized  model  can  now  be  generated  from  the  following  seven  flight 
parameters. 


M    Q    a    0    p    q    r 


(5.12) 


5J The  Flight  Parameter  Generator 

A  series  of  flight  conditions  are  made  and  used  to  generate  many  linear  models. 
Feedback  gains  are  generated  for  each  condition.  The  first  of  the  series  is  called 
the  nominal  flight  condition.  The  values  of  the  parameters  for  the  nominal  flight 
condition  are 

M  =  2,  Q  =  1250  psf,  a  =  8°,0  =  0°, 
p  =  07s,  ?=  107s,  r  =  07s. 

Next,  one  of  the  parameters  is  allowed  to  vary  while  the  other  six  are  held  constant. 
This  process  is  repeated  six  times  so  that  each  parameter  can  be  tested.  Table  5.1 
shows  the  starting  point,  and  the  minimum  and  maximum  values  of  each  changing 
parameter.  Parameters  with  nonzero  starting  points  begin  at  the  starting  point, 
increment  to  the  maximum  value,  return  to  the  starting  point,  and  then  decrement 
to  the  minimum  value.  Due  to  symmetry,  the  remaining  variables  have  starting  points 
at  zero  and  increment  to  their  maximum  values.  M  and  Q  sweep  through  a  narrow 
range  because  of  restrictions  of  the  atmospheric  tables. 

5J — Initializing  the  Iterative  Lvapunov  Design  Method 
The  iterative  Lyapunov  design  method  requires  an  initial  guess,  A'?0  and  Kpr0  and 
two  positive  definite  matrices,  Pq  and  Ppr.   The  initializer  supplies  these  values  and 
will  be  discussed  in  this  section. 
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Table 


Initial  point 

Minimum  value 

Maximum  value 

M 

2 

1 

2.6 

Q 

1250  psf 

700  psf 

5000  psf 

a 

8° 

-8° 

20° 

0 

0° 

0° 

10° 

P 

0°/s 

0°/s 

50075 

<? 

!0°/s 

-107s 

2007s 

r 

0°/s 

0°/s 

2007s 

The  initial  feedback  gains  are  found  by  using  a  pole  placement  algorithm.  At  the 
nominal  flight  condition,  the  linear  models  are  given  by 


A      — 

■n-pr  — 


Aq  — 


-0.459 

-2255.5 

73.0 


-1.1345    0.9996 
-261.4732    0.6209 


,-S,  = 


-0.1463 
-123.1091 


0.140  -.100 
-2.41  .066 
-0.181    -.648 


1  -Dpr   — 


0.018 

-1173.5 

2.01 


0.117 
-1335.6 

-114.4 


(5.13) 


(5.14) 


The  desired  eigenvalues  for  the  closed  loop  pitch  dynamics  have  been  chosen  to  be 
-40±jl0.  For  the  closed  loop  roll-yaw  model  the  desired  pole  locations  are  -20±j5. 
-80.  The  resulting  feedback  gains  are 


Kq  = 


-11.1633    -.6233 


,  J\pr  — 


-1.4166    -.0758       .3758 
2.9337       .0086    -.3300 


(5.15) 


For  both  closed  loop  systems,  P  must  be  found  so  that  xTPx  is  a  Lyapunov 
function.  The  following  problem  is  stated. 

Given  a  stable  linear  system  x  =  Ax,  find  a  positive  definite  function, 
V  =  xTPx  so  that  V  =  xT(PA  +  ATP)x  is  a  negative  definite  function. 

Let  A  be  put  into  Jordan  canonical  form. 


A  =  SJS 


-1 


(5.16) 
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where  S  is  an  invertible  matrix.  The  diagonals  of  J  are  the  real  parts  of  the  eigen- 
values of  A,  and  imaginary  parts  of  the  eigenvalues  lie  in  skew  symmetric  locations 
off  of  the  diagonals.  For  example,  if  the  eigenvalues  of  A  are  a  ±  jb,  c,  then 


J  = 


a 

-b    0  " 

b 

a    0 

0 

0     c_ 

(5.17) 


Let 


«*    —  Jdiaq  T  Js 


'skew 


(5.18) 

where  Jdmg  is  the  symmetric  part  of  J  and  Jskew  is  the  skew  symmetric  part  of  J.  In 
the  above  example 


** diag  — 


a 

0 

0  " 

0 

a 

0 

0 

0 

c 

*J skp.iu    — 


0 

-b    0  ' 

b 

0    0 

0 

0    0 

(5.19) 


Making  the  following  transformation  on  the  syst 


em  x  =  Ax,  let 


Then 


and 


x  =  Sz 


Si  =  ASz 


(5.20) 


(5.21) 


z  =  S  x ASz  =  Jz.  (5.22) 

Let  Pz  =  -Jdiag.  We  wish  to  check  the  velocities  of  the  system  z  =  Jz  on  the  ellipsoid 


zTPzz  =  -zTJdiagz  =  1 


(5.23) 


The  normal  to  the  ellipsoid  at  z  is  -JdiagZ.  The  projection  of  z  =  Jz  onto  the  normal 
is  JdiagZ.  Here  the  normal  velocity  component  has  the  same  magnitude  but  opposite 
direction  to  the  normal  of  the  ellipsoid.  If  the  system  has  no  complex  eigenvalues 
then  the  velocities  are  orthogonal  to  the  ellipsoid  everywhere.   For  this  reason,  the 
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choice  of  Pz  -  -Jdiag  is  the  best  choice  for  a  positive  definite  function  for  the  system 
z  =  Jz. 


V{z)  =  -zTJdiagz 
Making  the  following  transformation  into  the  x  coordinate  system  gives 


z  =  S~h 


(5.24) 


(5.25) 


which  implies 

V(x)  =  -xT[S-^jdmgS-^  (5.26) 

Our  choice  of  i3  is 

P  =  -[S-1}TJd*agS-K  (5.2T) 

For  the  nominal  flight  condition,  the  Jordan  canonical  form  of  Aq  -  BaKa  and 
Apr  —  BpTKpr  is  found  and  from  (5.27) 


'qJ-iq 


26772     609 
609    14.9 


p    — 

x  pr  — 


6792       10.7    -316.4 

10.7       4.02      -.509 

-316.4    -.509         15.7 


(5.28) 


The  eigenvalues  of  \(PqAq  +  ATqPq  -  PqBqKq  -  KjBjPq)  at  the  nominal  point 


are 


-1.0715  x  106,    -40.002 
Likewise,  for  the  closed  loop  roll-yaw  system  the  eigenvalues 


-1.3614  x  10s, -320.00, -20.003 


(5.29) 


are 


;s.30) 


<L4 — The  Iterative  Lvapunov  Design  Method 

The  iterative  Lyapunov  design  method  generates  feedback  gains  so  that  xTP  x 

and  xTPprx  are  Lyapunov  functions  for  each  closed  loop  system.    The  algorithm 

requires  the  initial  guesses  Kq0  and  A^o,  for  the  first  point,  and  the  positive  definite 

matrices  Pq  and  Ppr.    As  a  given  flight  condition  changes,  the  feedback  gains  from 
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the  previous  point  become  the  initial  guess  for  the  present  point.  The  result  of 
this  algorithm  is  a  series  of  gains;  one  set  for  each  flight  condition  generated.  The 
next  portion  of  this  section  discusses  some  modifications  made  to  the  algorithm  from 
section  2.4.  The  material  which  follows  describes  how  the  design  constraints  are 
selected. 

The  design  method  of  section  2.4  finds  a  K  so  that  the  eigenvalues  of  \{PA+ATP) 
fall  between  cf  and  cf  where  A  is  the  closed  loop  system.  The  design  algorithm 
used  in  Figure  5.1  has  been  modified  so  that  the  resulting  K  satisfies  the  following 
requirements. 


fii     <    Ax 

0.2       <       ^2 


\[P(A-  BK)  +  {A-  BK)T P] 
\[P(A-  BK)  +  (A-  BK)7 'P] 


1  I         <5-3r 


Qn  <  A„  [\[P(A  -  BK)  +  (A-  BK)TP] 
where  Aj  is  the  smallest  eigenvalue,  A2  is  the  second  smallest  eigenvalue,  and  so  on. 
The  values  of  the  c's  are  supplied  by  the  designer.  This  modification  has  been  made 
with  the  belief  if  more  design  constraints  are  made,  then  the  performance  will  vary 
less  with  changes  in  the  flight  conditions.  The  modified  algorithm  is  as  follows. 

1)  Choose  P,  clf,  clf,  ...  ,  cnf,  and  cn/.  The  selection  of  clf  and  cnf  must 
obey  Theorem  2.  However,  Theorem  2  does  not  guarantee  the  nonempti- 
ness  of  K,. 

2)  Compute  c10,c10,  ...  ,cn0,cn0  so  that  the  initial  guess,  KQ  €  £0  and 
rCf  C  JCq. 

3)  Let  i  =  0 

4)  Let  %  =  i  +  1 

5)  Compute  K%  so  that  Kt  e  £,-_i  and  Kx  ft  dfCt^. 

6)  Compute  cl2,  cu,  ...  cm,  cm  so  that  K{  £  did  and  Kj  C  Kt. 
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7)  Repeat  steps  4)  5)  and  6)  until  one  of  two  events  occur. 
L  £w  =  fii/,  cu  =  cif,  ...  ,  cnt  =  cnf  and  cni  =  cn/. 
2-  (fin  -Ci.-.j)  ...  (cni  -  cm_j)  become  very  small. 

Event  1.  implies  that  the  final  answer  has  been  found.  Event  2.  implies  that  JCf 
is  empty  and  Kt  satisfies  the  constraints  corresponding  to  cU:  clh  ...  ,  cm,  cm.  Steps 
2),  5),  and  6)  are  accomplished  in  the  following  way. 

Steps  2)  and  6) 


Let 


Then, 


el 


£-n        — 


Ai  [\{PA  +  ATP  -  PBK%  -  KJBTP) 
Xn  [\{PA  +  ATP  -  PBKi  -  KfBTP) 


(5.32) 


cu    = 


min[cx^  ei] 
max[cif,ei] 


=     min\c, 


nfi  enJ 


c-ni     =    max[cnf,en] 


(5.33) 


Step  5) 


Let  K  =  Ki 
For  j  =  1  to  m 

For  k  =  1  to  n 

P^-iPBQ^  +  Qj^P) 

P0  =  p(A  -  BKoj,k)  +  (A-  BKQhk)TP  -  2d 

Kvj,k  =  K  -  kj<kQjzk 
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*-2n- 


Solve  det\kPt  +  P0]  =  0  for  in  terms  of  k  for  c  =  c,,,  ...     L 
Reject  all  values  of  k  which  do  not  meet  the  constraints  from 
steps  2)  or  6).  The  result  is  2n  intervals  whose  lower  bounds  are 
designated  by  k^,  ...  ,  k2n  and  whose  upper  bounds  is  named  kx, 
Find  the  intersection  of  these  intervals  by  evaluating 
k  =  maa:[k1,...,k2n] 
k  =  mm[ki,...,k2„]. 
Find  the  midpoint  of  the  interval  by  computing 

kj,*  =  f(k  +  k). 
Let  K  =  kjikQj<k  +  KQjtk 
Next  k 
Next  j 
Ki+1  =  K 


Remarks: 

Let  the  feedback  sets  £l5  %u  £2,  JC2,  ...  ,  £n,  and  Kn  be  defined  respectively  as  the 
set  of  all  K  so  the  constraints  (5.31)  are  met.  Theorem  2  provides  conditions  for  the 
nonemptiness  of  £:  and  Xn.  But  conditions  for  the  nonemptiness  of  the  remaining 
sets  are  still  unknown.  Also,  Kt,  ...  JCn  are  not  convex  in  general.  These  are  the 
limitations  of  using  the  modified  design  algorithm. 

We  now  turn  to  selecting  constraining  values  for  the  eigenvalues  oi\[P{A-BK)  + 
{A  -  BK)TP].  It  is  necessary  to  evaluate  the  uncontrollable  normal  velocity  compo- 
nents for  each  flight  condition  that  will  be  generated  in  Figure  5.1.  Tables  5.2  and  5.3 
show  for  both  models  the  minimum  and  maximum  uncontrollable  normal  velocities 
for  each  changing  flight  parameter. 
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Table  5.2.  Uncontrollable  normal  velocity  components  for  the  pitch  model 


changing  parameter 


M 


Q 


a 


0 


P 


minimum  uncontrollable  velocity 


-43.1456 


-45.0144 


-43.3856 


-42.6047 


-42.6047 


-42.6047 


-42.6047 


maximum  uncontrollable  velocity 


-41.7253 


-42.2289 


-42.2362 


-42.6047 


-42.6045 


-42.6047 


-42.6044 


Table  5.3.  Uncontrollable  normal  velocity  components  for  the  roll-yaw  model 


changing  parameter 


M 


Q 


a 


(3 


minimum  uncontrollable  velocity 


-21.3937 


-22.3016 


-21.7136 


-21.3180 


-21.3181 


-21.3184 


-21.3180 


maximum  uncontrollable  velocity 


-20.9751 


-21.1574 


-21.1140 


-21.3180 


-21.3180 


-21.3180 


-21.3180 


For  the  pitch  model,  the  uncontrollable  normal  velocity  components  range  from 
-45.0144  to  -41.7253.  Theorem  2  requires  that 


-41.7253     < 


J?2 


s9l 


<     -45.0144 


(5.34) 


In  addition,  from  (5.29),  we  want 


cql    <    -1.0715  x  106    <    cql 
cq2    <         -40.002         <    cq2 


(5.35) 


From  this  the  constraining  values  for  the  pitch  model  have  been  chosen  to  be 


c?1  =  -1.2  x  106,    qc1  =  -lx  106, 


(5.36) 


«93  =•  -45,  cq2  =  -35 

Likewise  for  the  roll-yaw  model,  when  looking  at  Table  5.3,  Theorem  2  requires  that 


-20.9751     <         cpr3 
cPri  <    -22.3016 


(5.37' 
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The  eigenvalues  in  (5.30)  suggest  the  following. 

cprl     <     -1.3614  x  10s  <  cprl 

cpr2    <    -320.00    <  cpr2  (5.38) 

cpr3    <  -20.003    <  cpT3 

The  following  constraining  values  have  been  chosen  for  the  roll-yaw  model. 

cPri  =  -150000,    cprl  =  -110000 
cpr2  =  -340,  cpr2  =  -300  (5.39) 

£pr3  ~  ~^^i  Cpr3  =  — 17 

iL5 Formulation  of  a  State  Tracker 

The  autopilot  of  the  EMRAAT  missile  will  be  a  state  tracker.  That  is,  we  want  to 
be  able  to  change  the  location  of  the  equilibrium  point  in  order  to  control  the  values 
of  some  of  the  states.  The  following  shows  how  this  will  be  accomplished. 

Given  the  linear  system 

x  =  Ax  +  Bu  (5.40) 

y  =  C*i  (5.41) 

we  would  like  to  find  a  control  law 

u  =  -Kx  +  Krefv  (5.42) 

so  that  y,  the  output,  tracks  v  ,the  reference  input,  asymptotically.  We  require  y  =  v 
when  x  =  0.  When  x  =  0,  then 

0  =  ^x  -  BKx  +  BKrefv.  (5.43) 

Since  K  is  chosen  so  that  the  system  is  stable,  then  A  -  BK  is  invertible  and 

x  =  ~(A-BK)-lBKrefv  (5.44) 


Also 


V  = 


=  y  =  Cx  =  -C(A  -  BK)-lBKre;v  (5.45) 


72 


Because  (5.45)  is  true  for  any  v,  then 


I=-C{A-BK)-lBK, 


ref 


(5.46) 


Controllability  of  the  system  implies  that  C(A- BK)^B  is  invertible.  Solving  (5.46) 
for  Kref  gives 

Kref  =  -[C{A-BK)-lB}-'  (5.47) 

The  EMRAAT  missile  has  three  inputs  and  therefore  only  three  states  can  be 
tracked.  Controlling  a  in  the  pitch  model  and  fi  and  p  in  the  roll  yaw  model  is 
desirable.  For  the  pitch  model  y  =  a,  implying  that 


and,  therefore, 


For  the  roll-yaw  model 


Cg  =  [l    0] 


#«/,  =  ~{Cq(Aq  -  B.K^B,]-1 


y  = 


0 

P 


Opr  — 


1    0    0 
0    1    0 


(5.48) 


(5.49) 


(5.50) 


and,  thus, 

Krefpr   =  -[CpT(Apr  -  BprKp^Bpr]-1  (5.51) 

Krefq  and  Krefpr  are  computed  for  each  flight  condition  and  then  compared  along 
with  Kg  and  Kpr  to  the  flight  conditions. 

M — Comparing  Gains  with  Flight  Parameters 
A  series  of  gains  have  been  generated  as  a  function  of  different  flight  conditions. 
Each  flight  parameter  has  been  swept  through  a  range  of  points  while  the  remaining 
six  have  been  held  constant.  In  order  to  compare  different  gains  on  the  same  input 
it  has  been  decided  to  use  the  products  of  gains  and  their  corresponding  terms  at 
typical  values.    For  example  a  typical  value  of  a  is  8°.     So  we  set  a0  equal  to  8° 
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Table  5.4.  Extreme  values  and  range  of  the  pitch  channel  control  terms  as  individual 
flight  parameters  vary    Gains  depend  mostly  on  M,  Q,  and  a. 


M 

Q 

a 

(3 

"0  1    "•"^ 

p 

q 

r 

control 

min 

min 

min 

min 

min 

min 

min 

term 

max 

max 

max 

max 

max 

max 

max 

diff 

diff 

diff 

diff 

diff 

diff 

diff 

kqlt^O 

-2.15 

-3.1 

-2.1 

-1.60 

-1.61 

-1.61 

-1.60 

-.72 

-.16 

-1.4 

-1.59 

-1.59 

-1.59 

-1.59 

1.43 

2.95 

.69 

.017 

.017 

.017 

.017 

kql2<]0 

-.1509 

-.21 

-.16 

-.118 

-.118 

-.118 

-.118 

-.0529 

-.031 

-.10 

-.117 

-.117 

-.117 

-.117 

.099 

.179 

.054 

.0005 

.0005 

.0005 

.0005 

kref  qll&cO 

-2.53 

-3.5 

-2.7 

-1.98 

-1.98 

-1.98 

-1.98 

-.916 

-.53 

-1.7 

-1.97 

-1.96 

-1.96 

-1.96 

1.61 

3.0 

.94 

.017 

.017 

.017 

.017 

and  look  at  kmaQ.  We  set  q0  and  Qc0  equal  to  10°/s  and  8°  respectively  so  that 
we  can  look  at  kql2q0  and  kref  qUacQ.  The  sum  of  these  three  terms  are  fed  into  to 
the  elevator.    For  the  terms  which  are  fed  into  the  remaining  inputs,  the  following 

assignments  are  made. 

A)    =       &o       =         2° 

(5.52) 


Po    =      Pco      =    100% 
r0     =    25% 

Figures  5.2a-g  show  kqlla0  plotted  against  all  seven  flight  parameters.  These  seven 
figures  show  that  kqlla0  changes  with  M,  Q,  and  a  but  remains  nearly  constant  when 
/?,  p,  q,  and  r  change.  Similar  figures  exist  for  the  remaining  eleven  gains  and  are 
summarized  in  Tables  5.4  and  5.5.  The  minimum  and  maximum  values  of  each  term 
is  listed  for  each  changing  flight  parameter  along  with  the  difference  between  the 
minimum  and  maximum  values.  Angles  are  expressed  in  radians.  From  this  table  it 
was  determined  that  all  gains  will  be  scheduled  against  M,  Q,  and  a. 
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Figure  5.2.  kqlla0  vs.  a)Af;  b)Q;  c)q;  d)/5;  e)  p:  f)  ?; 
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Table  5.5.   Extreme  values  and  range  of  the  roll-yaw  channel  control  terms  as  indi- 

Vidua]   fllVht   naramptprc  irarv     n.-*;^,,   A 1 4.1 n  r     s\  i 


M 

Q 

a 

0 

P 

^  1    "-"■"■'■* 

9 

r 

control 

min 

min 

min 

min 

min 

min 

min 

term 

max 

max 

max 

max 

max 

max 

max 

diff 

diff 

diff 

diff 

diff 

diff 

diff 

KpriiPo 

-.050 

-.11 

-.22 

-.050 

-.050 

.050 

.050 

.022 

.045 

-.013 

-.038 

-.036 

-.036 

-.036 

.072 

.157 

.21 

.011 

.014 

.014 

.013 

Kprl2P0 

-.17 

-.14 

-.29 

-.13 

-.13 

-.13 

-.13 

-.011 

-.061 

-.13 

-.13 

-.13 

-.13 

-.13 

.16 

.20 

.15 

.001 

.001 

.001 

.001 

^prl3r0 

.074 

.047 

.15 

.16 

.16 

.16 

.16 

.18 

.28 

.25 

.16 

.17 

.17 

.17 

.10 

.23 

.099 

.0006 

.001 

.001 

.0009 

fc  pr  21 PO 

.026 

.0062 

.065 

.093 

.092 

.092 

.092 

.11 

.167 

.15 

.10 

.101 

.101 

.101 

.080 

.16 

.084 

.008 

.009 

.009 

.009 

kpr22P0 

-.054 

-.18 

.0015 

.014 

.014 

.014 

.014 

.015 

.015 

.152 

.014 

.015 

.015 

.015 

.068 

.20 

.15 

.0007 

.0008 

.0005 

.0007 

Kpr2Sr0 

-.151 

-.25 

-.24 

-.15 

-.15 

-.15 

-.15 

-.067 

-.034 

-.11 

-.14 

-.14 

-.14 

.14 

.084 

.22 

.12 

.001 

.002 

.001 

.002 

Kref  prllPcO 

-.14 

-.21 

-.16 

-.14 

-.14 

-.14 

-.15 

-.064 

-.050 

-.13 

-.13 

-.13 

-.13 

-.13 

.08 

.16 

.028 

.011 

.014 

.013 

.013 

Kref  prl2PcO 

-.071 

-.046 

-.39 

-.040 

-.040 

-.045 

-.040 

.029 

.15 

.056 

-.039 

-.040 

.013 

-.038 

.  .10 

.19 

.45 

.001 

.001 

.057 

.002 

"-re/  pr2lPcO 

.055 

.030 

.11 

.12 

.12 

.12 

.12 

.126 

.19 

.17 

.13 

.13 

.13 

.13 

.071 

.16 

.058 

.008 

.009 

.010 

.009 

Kref  pt22PcO 

-.010 

-.33 

-.17 

-.071 

-.071 

-.12 

-.071 

-.070 

-.070 

.24 

-.070 

-.070 

-.066 

-.070 

.027 

.26 

.41 

.001 

.001 

.051  1 

.002 

CHAPTER  6 
COMPUTING  LOOK-UP  TABLES 

In  the  last  chapter  we  showed  that  the  gains  depend  mostly  on  a,  M,  and  Q. 
This  chapter  describes  the  process  of  generating  a  look-up  table  of  gains  verses  M, 
Q,  and  a.  First  a  grid  of  points  is  formulated.  Design  constraints  are  then  formulated 
based  on  the  knowledge  of  uncontrollable  velocity  components  of  the  linear  models. 
Finally,  the  gains  are  computed. 

6T Determining  a  Grid  of  Points 

A  two  dimensional  grid  of  points  for  M  and  Q  has  been  made  and  used  for  each 
entry  of  a  in  the  table  for  both  the  pitch  channel  and  the  roll-yaw  channel.  Q  sweeps 
through  a  wide  range  of  values  starting  with  100  psf  and  ending  at  15,000  psf.  The 
values  of  M  were  chosen  so  that  each  entry  of  M  and  Q  lie  in  the  atmospheric  tables 
used  to  compute  p  and  V.  As  a  result,  the  grid  points  are  not  rectangular.  All 
entries  are  restricted  to  values  between  M  =  .6  and  M  =  3.5  and  must  correspond 
to  altitudes  between  sea  level  and  50,000  ft. 

Table  6.1  shows  the  values  of  a  used  in  the  look-up  table  for  both  channels.  M 
and  Q  sweep  through  all  values  of  the  grid  previously  mentioned  for  each  value  of  a 
in  Table  6.1.  The  spacing  between  the  grid  points  was  determined  in  a  trial  and  error 
process.  During  the  iterative  procedure  for  computing  feedback  gains,  the  initial 
guess  for  each  point  came  from  the  result  of  an  adjacent  grid  point.  The  closer  the 
spacing  between  adjacent  points,  the  fewer  iterations  were  needed  to  find  the  next 
feedback  gain.  Numerical  problems  as  described  in  section  2.4  and  shown  in  figure 
2.10  were  encountered.  When  this  happened  some  of  the  constraints  were  relaxed 
so  interior  points  in  the  feedback  set  could  be  found.    Later  these  constraints  were 
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Table  6.1.  Values  of  a  in  the  pitch  controller  look-up  tabk 


Pitch 


-3.0 


1.0 


4.0 


8.0 


12.0 


16.0 


20.0 


a 


Roll-Yaw 


-1.0 


1.0 


2.5 


4.0 


8.0 


12.0 


16.0 


20.0 


made  more  restrictive  and  returned  to  their  original  assignments.  This  became  a 
tedious  process  for  some  parts  of  the  grid.  When  the  number  of  iterations  exceeded 
100  it  was  decided  to  add  more  grid  points  so  that  the  desired  feedback  gains  could 
be  found  in  fewer  iterations. 

6J2 Formulation  of  the  Design  Constraints 

Before  using  the  iterative  Lyapunov  design  algorithm,  the  uncontrollable  normal 
velocity  components  for  the  entire  M-Q-a  grid  must  be  determined.  This  information 
is  needed  to  formulate  the  design  constraints.  This  process  is  described  by  the  block 
diagram  in  Figure  6.1.  Table  6.2  presents  the  minimum  and  maximum  uncontrollable 
normal  velocity  components  of  the  pitch  model  throughout  the  M  -  Q  grid  for  each 
value  of  a.  Here  P  =  Pq,  the  matrix  computed  during  the  initializing  procedure 
of  the  last  chapter.  Table  6.3  gives  the  same  result  for  the  roll-yaw  model  where 
P  =  PpT.  The  extreme  values  of  uncontrollable  normal  velocities  for  the  pitch  model 
are  -49.5985  and  -37.3480.  Also  the  eigenvalues  of 


\[Pq(Aq-BqKq)  +  (Aq-BqKqfPq] 


(6.i; 


pr 
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Figure  6.1.  Formulation  of  the  Design  Constraints. 
Table  6.2.  Uncontrollable  normal  velocity  components  for  the  pitch  model. 


r 


I — 7-j r-n — --  ~ — ™  ^uif  v^in,iL.y  cumponents  ior  tne  pitch  model. 

|  tt(degrees)  ||  minimum  uncontrollable  velocity  j  maximum  uncontrollable  velocity 

-37.4260 


1 


12 


16 


20 


-47.1539 


-47.0806 


-47.1611 


-47.9797 


-48.8796 


-49.5915 


-37.3873 


-37.3480 


-37.4247 


-37.3906 


-49.5985 


-37.6144 


-37.4355 


Table  6.3    Uncontrollable  normal  velocity  components  for  the  roll-yaw  model 


a( degrees) 


12 
16 
20 


minimum  uncontrollable  velocity 


-23.3399 


-23.3774 


maximum  uncontrollable  velocity 


-19.6307 


-23.4445 


-23.4971 


-24.5852 


-25.3163 
-25.6480 
-26.4797 


-19.6409 


-19.6883 


-19.7450 


-17.2696 


-19.9600 
-19.8625 
-20.3071 
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at  the  nominal  flight  condition  from  the  previous  chapter  were  found  to  be  -1.0715  x 
106  and  -40.002.  Theorem  2  requires  that 


-37.3480     <  cq2 

fill  <    -49.5985  (6-2) 


But  we  also  want 


%     <     -1.0715  x  106     <    cql 

Qq2     <  -40.002  <    cq2  (6-3) 

since  we  desire  the  eigenvalues  of  (6.1)  to  be  close  to  those  at  the  nominal  flight 

condition.  From  this,  the  constraining  values  were  chosen  to  be 

cql  =  -1.2  x  106,  cql  =  -106, 

cq2  =  -45,  cq2  =  -35 

Likewise,  for  the  roll-yaw  model,  the  uncontrollable  normal  velocity  components  range 
from  -26.4797  to  -17.2696.  From  Theorem  2  we  must  have 

-17.2696    <         cpr3 

cpri         <    -26.4797  (6-4) 

With  the  eigenvalues  of 


^[PPr(Apr  -  BprKpr)  +  (APT  -  BpTKpr)TPpr] 

at  the  nominal  flight  condition  being  -1.3614  x  105,   -320.00,  and  -20.003,  the 
following  is  desirable. 

cprl     <    -1.3614  x  105  <  cpTl 

cpr2    <         -320.00    <  cpr2  (6.5) 

Qpr3    <         -20.003    <  cpr3 

With  this  in  mind,  the  following  selections  were  made. 

cprl  =  -200000,  cprl  =  -90000, 

cpr2  =  -400,  cpr2  =  -250, 

£Pr3  =  —25,  cpr3  =  -15 
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6^3 Generating  the  Look-Up  Table 

With  the  design  constraints  set,  feedback  gains  and  feedforward  gains  are  gen- 
erated for  each  grid  point.  Figure  6.2  gives  a  block  diagram  of  the  system  used  to 
accomplish  this.  For  the  pitch  model,  the  initial  guess  comes  from  one  of  the  gains 
that  was  generated  in  the  previous  chapter  when  determining  the  dependence  be- 
tween gains  and  flight  parameters.  The  operating  point  from  which  this  initial  guess 
originates  is 

M  =  1,  Q  =  1250psf,  a  =  8°  (6.6) 

and  is  one  of  the  extreme  values  listed  in  Table  5.1.  The  result  of  the  first  point  is 
used  to  start  adjacent  points  which,  in  turn,  start  new  adjacent  points  until  gains 
have  been  computed  for  the  entire  grid.  The  look-up  table  for  the  roll-yaw  model  is 
made  in  the  same  way. 

Numerical  problems  were  encountered  in  parts  of  the  look-up  table  for  the  roll- 
yaw  model.  They  were  similar  to  the  problems  that  were  predicted  in  step  5)  of  the 
iterative  design  method  given  in  Chapter  2.  To  overcome  these  difficulties,  some  of  the 
constraining  values  were  relaxed  for  a  number  of  iterations  and  were  later  returned  to 
their  original  assignments  in  the  algorithm.  Eventually,  the  desired  feedback  matrix 
was  found. 

Figure  6.3  shows  km  verses  mach  number  and  dynamic  pressure  when  a  =  8°.  As 
this  figure  would  indicate,  the  gains  generated  from  the  iterative  Lyapunov  method 
are  smooth  with  respect  to  the  dependent  flight  parameters.  This  fact  gives  hope 
that  the  gain  scheduling  scheme  will  be  easy. 
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Figure  6.2.  Formulation  of  the  look-up  table. 
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Figure  6.3.  kqll  verses  M  and  Q  when  a  =  8° 


CHAPTER  7 
GAIN  SCHEDULING 

The  result  of  the  previous  chapter  is  thirteen  tables  of  gains  in  terms  of  mach 
number,  dynamic  pressure,  and  angle  of  attack.  This  chapter  discusses  the  process 
of  curve  fitting  used  to  implement  the  look-up  table.  The  results  of  a  test  of  this 
scheme  will  follow. 

7.1     Curve  Fitting 

It  was  decided  to  use  a  combination  of  polynomial  fitting  and  interpolation  to 
implement  the  autopilot.  Third  order  polynomials  were  fit  to  the  tables  as  a  function 
of  mach  number.    However,  low  order  polynomials  could  not  achieve  close  fits  as  a 
function  of  angle  of  attack  or  dynamic  pressure,  so  linear  interpolation  was  used  for 
these  two  variables.   Figure  7.1a  shows  kqll  as  a  function  of  M  for  various  constant 
values  of  Q  at  a  =  8°.    An  example  of  polynomial  fitting  of  one  of  these  curves  is 
shown  in  Figure  7.1b  where  Q  =  150psf.    A  polynomial  has  been  made  for  every 
a-Q  pair  and  the  polynomial  coefficients  are  interpolated  as  a  function  of  these  two 
variables.    The  tables  of  the  three  pitch  controller  gains  each  have  7  entries  for  a 
and  22  entries  for  Q.   Since  each  polynomial  has  4  coefficients  the  total  number  of 
coefficients  for  each  gain  is  7  x  22  x  4  =  616.    Similarly  the  tables  for  the  roll-yaw 
controller  have  8  entries  for  a  and  the  same  number  of  entries  for  Q.  Each  of  the  ten 
gains  are  then  scheduled  using  8  x  22  x  4  =  704  polynomial  coefficients. 

7.2     Testing  the  Fit 

Figure  7.2  gives  the  system  for  testing  the  polynomial  and  interpolation  routines. 
Gains  were  generated  from  these  routines  at  locations  centered  between  the  original 
grid  points.  These  new  locations  were  found  by  taking  the  average  of  the  coordinates 

83 


84 


Figure  7.1  A  plot  of  the  kqll  verses  M  with  (a)  various  constant  values  of  Q  and 
a  -  8  ,  and  (b)  a  least  squares  third  order  polynomial  fit  for  the  plot  where  0  = 
150psf  and  a  =  8°. 

of  adjacent  grid  points.  Linear  models  were  also  made  at  each  test  point  where  /?,  p, 
?,  and  r  were  set  to  zero.  The  eigenvalues  of 


-z[Pq(Aq-BqKq)  +  (Aq-BqKq)TPq] 


(7.1) 


and 


^[PPAAr  -  BprKpr)  +  (Apr  -  BprKpT)TPpr]  (7.2) 

were  computed  to  see  if  they  remained  within  the  desired  limits.  Table  7.1  shows 
the  maximum  eigenvalue  of  (7.2)  for  all  of  the  test  points  at  a  =  1.75°.  These  values 
are  plotted  against  their  indices  in  Figure  7.3.  Most  of  the  eigenvalues  of  Table 
7.1  lie  within  the  desired  limits  of  -25  and  -15.  The  eigenvalue  in  seventeenth  row, 
second  column,  however,  is  -3.96,  the  worst  value  found  out  of  all  of  the  test  points. 
Although  the  deviation  is  high,  this  value  is  still  negative  indicating  stability  for  the 
closed  loop  system.  For  the  pitch  controller,  the  actual  limiting  values  of  A:  and  A2 
at  all  the  test  points  are 


-1.37  x  106    <    A?j     <    -9.85  x  10s 
-47.7         <    \q2    <         -38.1 


7.3: 


85 


Flight 

Parameter 

Generator 


a 


M 


Q 


Po- 


Linear 
Model 


V 


Atmospheric 
Tables 


B„ 


A 


pr 


B,n 


Polynomial  and 

Interpolation 

Routines 


Ka 


K, 


pr 


Kr 


e/» 


K 


reft 


£5. 


MVa 


A  (VT. 


Pr  ; 


P   P 


\ 


A's 


Comparison 
with  Design 
Constraints 


Figure  7.2.  A  test  of  the  curve  fitting  routines  used  to  implement  the  autopilot. 
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Table  7.1.  The  maximum  eigenvalues  of  (6.1)  for  all  test  points  at  a  =  1.75" 
Q/M  indices 


6 


9 


10 


11 


12 


13 


14 


15 


16 


17 


18 


19 


20 


21 


2 


-20.79 


-20.68 


-20.84 


-20.59 


-20.30 


-20.41 


-20.82 


-20.69 


-20.85 


-20.61 


-20.22 


-20.84 


-20.70 


-20.84 


-20.49 


-19.97 


-19.10 


-18.62 


-19.90 


-20.23 


-20.57 


-20.77 


-19.83 


-20.45 


-22.42 


-21.82 


-21.85 


-18.88 


-20.85 


-21.65 


-20.62 


-17.70 


-19.79 


-20.97 


-20.88 


-20.38 


-17.82 


-18.51 


-22.23 


-22.13 


-19.45 


-3.96 


-19.74 


-21.06 


-21.93 


-22.15 


-19.93 


-19.77 


-20.87 


-20.72 


-20.83 


-20.30 


-20.88 


-20.74 


6 


-20.89 


-20.76 


-20.81 


-20.18 


-19.71 


-17.86 


-20.67 


-19.33 


-11.67 


-12.40 


-17.66 


-21.91 


-20.16 


-18.92 


-20.37 


-14.48 


-16.89 


-18.26 


-20.97 


-19.96 


-19.88 


-20.03 


-16.16 


-20.79 


-20.33 


-19.88 


7 


-20.90 


-20.76 


-20.78 


-20.34 


-20.01 


-20.90 


-20.76 


-20.79 


-20.21 


-20.20 


-20.42 


-19.58 


-19.96 


-20.60 


-20.57 


-18.83 


-15.17 


-20.63 


-21.48 


-21.71 


-22.29 


-20.81 


-17.42 


-18.24 


-20.54 


-22.26 


-22.29 


-22.52 


-22.27 


-21.85 


-20.87 


-20.84 


-20.65 


-22.63 


-22.62 


-21.09 


-19.14 


-21.08 


-20.40 


-20.85 


-21.65 


-20.79 


-19.75 


-20.62 


-20.97 


-20.21 


-20.93 


-19.72 


-20.28 


-22.17 


-22.21 


-22.35 


-22.58 


-22.69 


-22.93 


-21.23 


-20.92 


-22.13 


-22.65 


-22.98 


-21.14 


-22.13 


-23.19 


-22.61     -22.53 


-22.42 


-21.49 


-20.94 
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Figure  7.3.  A  plot  of  the  eigenvalues  from  Table  6.1. 


87 

Likewise,  for  the  roll-yaw  controller,  the  limiting  values  at  all  the  test  points  are 

-2.03  x  105    <    Xprl     <    -9.28  x  104 
-2392.7        <    Xpr2    <        -182.6  (74) 

-26.27        <    Apr3    <         -3.96 

Some  of  these  values  differ  significantly  from  the  desired  constraints;  however,  since 
these  values  are  still  negative,  these  deviations  are  acceptable  and  indicate  that  the 
closed  loop  system  will  be  stable.  It  should  also  be  noted  that  most  eigenvalues 
remain  well  within  their  desired  constraints  as  shown  in  Table  7.1. 


CHAPTFR  8 
NONLINEAR  SIMULATIONS 

A  nonlinear  simulation  has  been  used  to  test  the  proposed  autopilot  for  the  EM- 
RAAT  missile.  First  a  section  follows  giving  an  overview  of  the  nonlinear  simulation. 
A  test  module  is  then  made  to  generate  state  commands  in  order  to  evaluate  the 
autopilot's  tracking  ability.  Finally,  a  series  of  flight  scenarios  are  run  to  determine 
the  ability  the  missile  has  to  intercept  the  target. 

8.1     The  Nonlinear  Simulation 

Figure  8.1  shows  a  block  diagram  of  the  simulation  used  to  test  the  EMRAAT 
missile.  The  program  is  written  in  FORTRAN.  Initial  conditions  of  the  target  and 
missile  are  specified  by  the  user.  The  simulation  is  then  run  and  a  trajectory  of  both 
the  target  and  missile  results.  All  target  and  missile  variables  can  be  observed.  The 
target  is  programmed  to  fly  in  a  straight  line  until  the  range  between  the  target  and 
missile  falls  below  5,000  ft.  The  target  then  makes  a  9  g  turn  to  the  right.  The 
simulation  terminates  when  the  closing  velocity  becomes  positive. 

The  seeker  measures  the  line  of  sight  angles  and  the  range  rate  of  the  target. 
The  simulation  uses  exact  measurements  of  these  values  and  does  not  assume  any 
noise.  These  values  are  sent  to  the  guidance  law  which,  in  this  case,  implements 
proportional  navigation.  A  derivation  of  this  guidance  law  can  be  found  in  Bryson 
and  Ho  [19].  The  outputs  of  the  guidance  law  are  two  desired  accelerations,  ayc  and 
aZc.  The  BTT  logic  makes  the  conversion  from  the  acceleration  commands  to  the 
three  state  commands  ac,  0C,  and  pc.  Since  the  missile  can  achieve  a  much  higher 
acceleration  with  angle  of  attack  than  with  sideslip,  the  BTT  logic  uses  pc  to  rotate 
the  desired  accelerations  into  the  pitch  plane.  If  this  roll  maneuver  is  successful  then 
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Figure  8.1.  A  Block  Diagram  of  the  Nonlinear  Missile  Simulation. 
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a*  will  become  small  and  aZc  will  become  positive.    ac  and  &  are  computed  in  an 
attempt  to  match  aZc  and  aVc  respectively. 
The  autopilot  implements  the  control  law 

u  =  -Kx  +  Krefv  (-81) 

where  x  is  a  vector  containing  the  actual  states  and  v  contains  the  state  commands 
from  the  bank-to-turn  (BTT)  logic.  The  states  come  from  exact  measurements  in  the 
simulation.  If  this  autopilot  were  to  be  implemented  in  an  actual  missile,  the  states 
would  be  measured  using  an  inertial  platform.  The  gains  *  and  Kref  come  from  the 
gain  schedule  implemented  with  a  combination  of  polynomials  and  interpolation.  In 
this  simulation  there  is  no  delay  in  the  gain  schedule  and  K  and  KTef  are  produced 
instantaneously.  The  output  of  the  autopilot  is  the  control  surface  angles  Sp,  Sq,  and 
ST.  Linear  and  angular  accelerations  are  computed  by  the  missile  dynamics  module 
of  the  program.  The  simulation  uses  the  output  of  the  missile  dynamics  to  compute 
all  of  the  flight  variables  including  the  position  and  velocity  of  the  missile. 

8.2  A  Test  of  State  Tracking 
The  model  for  the  EMRAAT  missile  has  five  states  and  three  inputs.  The  autopi- 
lot is  designed  to  track  three  state  commands:  aa,  &,  and  pc.  Before  running  missile 
target  scenarios  it  was  decided  to  test  the  autopilot's  tracking  ability.  The  BTT  logic 
was  disconnected,  and  the  following  commands  were  applied  to  the  reference  inputs 
of  the  autopilot. 

10°     for  0  s  <  t  <  .5  s 

0°      for  .5  s  <  t  <  2.75  s 

10°    for  2.75  s  <  t  <  3.75  s  (8"2) 

0°      for  3.75  s  <  t 

0°    for  0  s  <  t  <  2  s 
0c  =  {   5°    for  2  s  <  t  <  2.5  s  (g  3) 

0°    for  2.5  s  <  t 
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Table  8.1.  The  rise  times  of  each  commanded  state. 


Altitude(ft) 

Mach 

tra(s) 

trg(s) 

*rP(s) 

20,000 

2.0 

.100 

.243 

.0284 

50,000 

3.0 

.108 

.316 

.0181 

Pc  =    1 


f  0%  for  0  s  <  t  <  1  s 

1007s  for  1  s  <  t  <  1.5  s 

0°/s  for  1.5  s  <  t  <  2.75  s 

100%  for  2.75  s  <  t  <  3.25  s  (8-4) 

-100%  for  3.25  s  <  t  <  3.75  s 

I  0%  for  3.75  s  <  t 

The  experiment  was  performed  once  at  an  initial  altitude  of  20,000  ft  with  an  initial 
mach  number  of  2.0  and  a  second  time  at  an  initial  altitude  of  50,000  ft  with  an  initial 
mach  number  of  3.0.  Figure  8.2  shows  the  results  of  the  first  test.  In  addition  to  the 
commanded  states,  figure  8.2  also  shows  the  remaining  two  states  q  and  r.  The  rise 
time  as  defined  as  the  time  needed  to  achieve  90%  of  the  desired  value  was  found 
for  each  commanded  state  and  is  shown  in  table  8.1.  It  should  be  noted  that  mach 
number  and  dynamic  pressure  do  not  change  significantly  during  this  test.  This  is  a 
test  of  accuracy  in  state  tracking  and  is  not  a  valid  test  of  gain  scheduling  in  terms 
of  M  and  Q.  However,  a  changes  very  rapidly  and  does  not  appear  to  hinder  the 
performance. 

Cross  coupling  is  evident  from  this  experiment.  The  step  in  the  roll  rate  of  figure 
8.2c  at  t  =  2  s  is  due  to  the  sideslip  command  in  figure  8.2b.  A  lesser  degree  of  cross 
coupling  occurs  during  the  roll  command  for  2.75  s  <  t  <  3.75  s  when  a  which  is  at 
10°  is  rotated  into  sideslip.  The  same  effects  are  present  in  the  second  experiment  at 
50,000  ft. 

Figure  8. 2d  and  8.2e  show  the  two  states  that  have  no  commands.  The  spikes  in 
the  pitch  rate  occur  when  ac  changes  value,  because  a  pitching  maneuver  is  required 
to  change  the  angle  of  attack.    The  two  spikes  in  the  yaw  rate  occur  for  the  same 
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Figure  8.2.  A  test  of  state  tracking  at  an  initial  altitude  of  20,000  ft  with  an  initial 
mach  number  of  2.0.  The  states  shown  are  (a)  a,  (b)  /?,  (c)  p,  (d)  q,  and  (e)  r. 
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reason  when  sideslip  changes.    The  two  steps  in  the  yaw  rate,  however,  are  due  to 
cross  coupling  with  the  roll  rate. 

The  rise  times  in  table  8.1  are  small  and  indicate  that  the  missile  will  perform  well 
in  flight  scenarios.  While  cross  coupling  effects  are  noticeable,  they  are  not  believed 
to  be  great  enough  to  significantly  hinder  the  performance  of  the  missile;  thus  it  was 
decided  to  run  the  missile  through  a  series  of  flight  scenarios  to  see  how  well  the 
missile  can  intercept  the  target. 

8J Simulation  of  Flight  Scenarios. 

A  series  of  flight  scenarios  has  been  run  to  test  the  autopilot.  Figure  8.3  shows 
the  trajectory  of  the  missile  and  target  for  a  flight  scenario  at  an  altitude  of  20,000 
ft.  The  miss  distance  is  .64  ft.  A  hit  is  considered  to  be  any  miss  distance  under  10 
ft.  Figure  8.4  shows  the  commanded  and  the  actual  y  and  z  accelerations.  The  devi- 
ations from  the  commanded  normal  accelerations  are  mostly  due  to  a  simplification 
used  in  implementing  the  BTT  logic.  Instead  of  using  two  aerodynamic  coefficients, 
proportionality  constants  were  assumed.  The  errors  are  not  large  enough  to  prevent 
the  interception  of  the  target.  The  state  tracking  of  the  reference  inputs  appears  to 
be  working  well  as  seen  in  figure  8.5  and  indicates  that  the  autopilot  is  performing 
well. 

The  trajectories  of  a  40,000  foot  altitude  scenario  is  shown  in  figure  8.6.  The  miss 
distance  is  0.05  ft.  A  similar  scenario  at  10,000  ft  is  shown  in  figure  8.7.  The  miss 
distance  is  2.8  ft. 

Figure  8.8  gives  a  case  were  a  miss  occurred.  The  scenario  took  place  at  50,000 
ft  and  the  target  was  missed  by  452  ft.  Figure  8.9a  and  b  shows  that  the  missile  was 
unable  to  achieve  the  desired  z  acceleration  and  angle  of  attack  after  1.3  seconds  into 
the  flight.  This  is  because  the  elevator  reached  its  -40°  limit  as  figure  8.9c  indicates. 
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Figure  8.3.  (a)  Top  view  and  (b)  side  view  of  missile  and  target  trajectories  of  a 
scenario  which  occurred  at  20,000  ft.  The  miss  distance  is  0.64  ft.  Initial  mach 
numbers  of  the  target  and  missile  are  2.5  and  .92  respectively. 
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a)  Commanded  z  Acceleration  (-)  and  Actual  z  Acceleration  (-) 
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Figure  8.4.    Commanded  and  actual  (a)  z  accelerations  and  (b)  y  accelerations  for 
the  scenario  in  figure  8.3. 
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a)  Commanded  Alpha  (-)  and  Actual  Alpha  (-) 
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Figure  8.5.   Commanded  and  actual  states,  (a)  a,  (b)  0,  and  (c)  p.  for  the  scenario 
in  figure  8.3. 
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Figure  8.6.  (a)  Top  view  and  (b)  side  view  of  missile  and  target  trajectories  of 
a  scenario  which  occurred  at  40,000ft.  The  miss  distance  is  .05  ft.  Initial  mach 
numbers  of  the  target  and  missile  are  3.0  and  .92  respectively. 
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Figure  8.7.  (a)  Top  view  and  (b)  side  view  of  missile  and  target  trajectories  of 
a  scenario  which  occurred  at  10,000ft.  The  miss  distance  is  2.8  ft.  Initial  mach 
numbers  of  the  target  and  missile  are  2.0  and  .92  respectively. 
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The  miss  occurred  due  to  the  reduced  elevator  effectiveness  from  the  low  air  density 
at  high  altitudes. 

Another  miss  is  shown  in  figure  8.10.  Figure  8.11  shows  that  the  desired  z  ac- 
celeration was  not  achieved  because  the  commanded  angle  of  attack  was  at  its  18° 
limit.  The  commanded  normal  acceleration  was  too  great  to  be  realized.  Again,  the 
autopilot  performed  as  well  as  it  could,  but  the  miss  occurred  because  of  physical 
limitations  of  the  missile.  The  results  of  these  simulations  indicate  that  the  autopilot 
does  its  job  and,  in  reasonable  scenarios,  allows  the  missile  to  intercept  the  target. 
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Figure  8.8.  (a)  Top  view  and  (b)  side  view  of  missile  and  target  trajectories  of  a 
scenario  which  occurred  at  50,000  ft.  The  missile  missed  the  target  by  452  ft.  Initial 
mach  numbers  of  the  target  and  missile  are  3.0  and  .92  respectively. 
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a)  Commanded  z  Acceleration  (-)  and  Acrual  z  Acceleration  (-) 
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Figure  8.9.  Commanded  and  actual  (a)  accelerations  in  the  z  direction  (b)  com- 
manded and  actual  angles  of  attack,  and  (c)  the  elevator  deflection  angle  for  the 
scenario  shown  in  figure  8.8.    Saturation  of  the  elevator  angle  was  the  cause  of  the 

miss. 
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i)  Top  View  of  Missile  (-)  and  Target  (-)  Trajectory 
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Figure  8.10.  (a)  Top  view  and  (b)  side  view  of  missile  and  target  trajectories  of  a 
scenario  which  occurred  at  35,000ft.  The  target  was  missed  767  ft.  Initial  mach 
numbers  of  the  target  and  missile  are  1.9  and  .92  respectively. 
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a)  Commanded  z  Acceleration  (-)  and  Actual  z  Acceleration  (-) 
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Figure  8.11.  Commanded  and  actual  accelerations  (a)  in  the  z  direction  (b)  com- 
manded and  actual  angles  of  attack,  and  (c)  elevator  deflection  angle  for  the  scenario 
shown  m  figure  8.10.  Saturation  of  the  a  command  resulted  in  the  miss 


CHAPTER  9 
CONCLUSION 

A  controller  design  algorithm  has  been  made  which  can  stabilize  time  varying 
linear  systems.  The  motivation  behind  this  work  was  the  need  to  achieve  local  sta- 
bility of  a  nonlinear  system  by  stabilizing  a  linearized  model  which  is  a  function  of 
time  varying  system  parameters.   This  work  began  by  studying  second  order  linear 
time  varying  systems  that  were  stable  when  frozen  at  any  given  time  but  were  un- 
stable when  allowed  to  vary  with  time.  Geometric  techniques  were  used  to  rederive 
Lyapunov's  linear  stability  theorem.  These  techniques  were  carried  further  to  give  a 
condition  for  the  existence  of  state  feedback  which  places  an  upper  and  lower  bound 
on  the  eigenvalues  of  the  derivative  of  a  given  Lyapunov  function.  This  led  to  a  design 
algorithm  which  makes  the  derivative  of  the  Lyapunov  function  negative  definite.  A 
time  varying  second  order  system  was  then  stabilized  when  ordinary  pole  placement 
techniques  failed.  The  algorithm  was  then  applied  to  a  linearized  model  of  the  EM- 
RAAT  missile  and  used  to  formulate  an  autopilot.    Feedback  gains  were  generated 
along  a  grid  of  flight  parameters.  The  gains  were  scheduled  against  angle  of  attack, 
mach  number,  and  dynamic  pressure.  The  resulting  autopilot  was  tested  using  a  non- 
linear simulation.  The  missile  was  able  to  intercept  the  target  in  reasonable  scenarios 
and  flight  conditions. 

The  significance  of  the  proposed  design  algorithm  is  that  no  assumption  is  made 
on  time  invariance  of  the  system.  However  as  stated  in  chapter  2,  this  algorithm 
can  only  be  used  to  stabilize  a  system  if  there  exists  a  constant  positive  definite  P 
such  that  Bl(t)[A(t)P^  +  P-1AT{t)]BL(t)  is  negative  definite  for  all  time.  If  this 
is  true  then  a  K(t)  can  be  found  so  that  xTPx  is  a  Lyapunov  function  for  the  time 
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varying  system.  Under  this  restriction,  the  algorithm  can  be  used  to  achieve  local 
stability  for  a  nonlinear  system  by  formulating  a  linearized  model  around  a  time 
varying  operating  point.  Using  this  model,  feedback  gains  are  then  computed  as  a 
function  of  parameters  which  characterize  the  operating  point.  There  is  no  restriction 
on  the  rate  of  change  of  the  system  parameters. 

In  stabilizing  an  time  varying  linear  system,  it  is  not  known  how  to  choose  a 
constant  positive  definite  P,  or  even  if  a  P  exists  so  Bl(t)[A(t)P~l  + P^AT(t)}Bj_(t) 
is  positive  definite.  Further  study  in  this  area  would  be  warranted. 


APPENDIX  A 
AERODYNAMIC  DATA  FOR  THE  EMRAAT  AIRFRAME 

This  appendix  gives  aerodynamic  data  for  the  extended  medium  range  air-to-air 
technology  (EMRAAT)  airframe.  The  aerodynamic  coefficients  are  given  in  tabular 
form.  All  coefficients  that  correspond  to  angle  rates  have  no  dimensions  and  must  be 
multiplied  by  £.  All  coefficients  that  correspond  to  angle  positions  are  given  in  per 
degrees  and  must  be  multiplied  by  i§2. 


Missile  Reference  Diameter 
Missile  Reference  Area 


d  =  .625  ft 
S  =  .3067  ft2 


Table  A.l.  Tabular  data  for  C 


N 


a/M 

.8 

1.2 

2.0 

2.7 

3.5 

-8° 

-4.816 

-5.472 

-5.104 

-4.502 

-3.952 

-4° 

-2.408 

-2.736 

-2.552 

-2.251 

-1.976 

0° 

0.000 

0.000 

0.000 

0.000 

0.000 

2° 

1.204 

1.368 

1.276 

1.126 

.988 

4° 

2.408 

2.736 

2.552 

2.251 

1.976 

8° 

5.264 

6.072 

5.644 

5.011 

4.456 

12° 

8.613 

10.068 

9.456 

8.458 

7.528 

16° 

12.467 

14.674 

13.898 

12.490 

11.092 

20° 

16.800 

19.730 

18.610 

16.732 

14.860 
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Table  A.2.  Tabular  data  for  C„ 


a/M 

.8 

1.2 

2.0 

2.7 

3.5 

-8° 

6.976 

12.336 

11.376 

8.445 

7.472 

-4° 

3.488 

6.168 

5.688 

4.222 

3.736 

0° 

0.000 

0.000 

0.000 

0.000 

0.000 

2° 

-1.746 

-3.084 

-2.844 

-2.111 

-1.868 

4° 

-3.488 

-6.168 

-5.688 

-4.222 

-3.736 

8° 

-7.762 

-11.238 

-10.662 

-8.174 

-6.908 

12° 

-11.802 

-15.530 

-15.074 

-11.439 

-9.972 

16° 

-14.880 

-19.330 

-19.282 

-14.603 

-13.624 

20° 



-16.030 

-25.700 

-29.560 

-24.218 

-25.930 

Table  A. 3.  Tabular  data  for  C 


N, 


6q 


a/M 

.8 

1.2 

2.0 

2.7 

3.5 

-8° 

.153 

.195 

.095 

.077 

.048 

-4° 

.181 

.215 

.106 

.073 

.055 

0° 

.197 

.219 

.107 

.069 

.056 

2° 

.199 

.212 

.104 

.068 

.055 

4° 

.198 

.207 

.101 

.067 

.051 

8° 

.199 

.195 

.091 

.056 

.046 

12° 

.199 

.175 

.076 

.047 

.044 

16° 

.188 

.161 

.064 

.044 

.047 

20° 

.184 

.150 

.056 

.047 

.055 
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Table  A.4.  Tabular  data  for  C„ 


•Sq 


a/M 

.8 

1.2 

2.0 

2.7 

3.5 

-8° 

-1.035 

-1.278 

-.5951 

-.5269 

-.3970 

-4° 

-1.200 

-1.44 

-.6862 

-.5490 

-.4377 

0° 

-1.320 

-1.4719 

-.7089 

-.5297 

-.4151 

2° 

-1.330 

-1.4875 

-.7008 

-.5269 

-.3925 

4° 

-1.330 

-1.456 

-.6862 

-.5297 

-.3849 

8° 

-1.330 

-1.325 

-.6276 

-.4772 

-.3396 

12° 

-1.330 

-1.216 

-.5626 

-.4200 

-.3245 

16° 

-1.286 

-1.114 

-.5100 

-.3641 

-.3547 

20° 

-1.232 

-1.0594 

-.4715 

-.3724 

-.4226 

Table  A. 5.  Tabular  data  for  Cn 


a/M 

.8 

1.2 

2.0 

2.7 

3.5 

-8° 

1.0702 

.8632 

.2175 

-.0923 

-.3292 

-4° 

1.2877 

1.2182 

.5090 

.1754 

-.0234 

0° 

1.3544 

1.3825 

.6372 

.3138 

.1846 

2° 

1.3544 

1.4414 

.6500 

.3446 

.2338 

4° 

1.3509 

1.4772 

.6372 

.3723 

.2677 

8° 

1.2561 

1.5172 

.491 

.3200 

.2862 

12° 

.9502 

1.4133 

.2385 

.2708 

.3169 

16° 

.6989 

1.1228 

.3808 

.4338 

.5046 

20° 

.5123 

.5179 

.6179 

.6554 

.7446 
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Table  A.  6.  Tabular  data  for  CYr 


a/M 

.8 

1.2 

2.0 

2.7 

3.5 

-8° 

-.2226 

-.2052 

-.2175 

-.1717 

-.1758 

-4° 

-.2870 

-.2668 

-.2305 

-.2195 

-.1911 

0° 

-.3096 

-.3099 

-.2546 

-.2478 

-.2096 

2° 

-.3210 

-.3235 

-.2660 

-.2579 

-.2299 

4° 

-.3200 

-.3363 

-.2717 

-.2629 

-.2478 

8° 

-.3113 

-.3614 

-.2851 

-.2623 

-.2812 

12° 

-.2918 

-.3774 

-.3238 

-.2925 

-.3178 

16° 

-.2445 

-.3447 

-.3849 

-.3698 

-.3783 

20° 

-.2515 

-.2730 

-.4800 

-.4440 

-.4631 

Table  A. 7.  Tabular  data  for  C, 


n6r 


a/M 

.8 

1.2 

2.0 

2.7 

3.5 

-8° 

-1.0842 

-.9703 

-.4520 

-.4236 

-.2831 

-4° 

-1.0743 

-.9505 

-.4600 

-.3709 

-.2739 

0° 

-1.0446 

-.9683 

-.4900 

-.3818 

-.3015 

2° 

-1.0356 

-.9802 

-.4900 

-.4000 

-.3162 

4° 

-1.0257 

-1.0050 

-.5080 

-.4091 

-.3272 

8° 

-1.0000 

-1.0535 

-.5200 

-.4436 

-.3640 

12° 

-.9505 

-1.0891 

-.5300 

-.5200 

-.4412 

16° 

-.8851 

-1.0931 

-.5760 

-.6400 

-.5441 

20° 

-.7921 

-1.0475 

-.6840 

-.7709 

-.6838 

110 


Table  A. 8.  Tabular  data  for  Q, 


a/M 

.8 

1.2 

2.0 

2.7 

3.5 

-8° 

-.11703 

-.11700 

-.08247 

-.06413 

-.04209 

-4° 

-.12446 

-.12390 

-.08446 

-.05906 

-.04209 

0° 

-.13189 

-.12870 

-.08705 

-.06123 

-.04420 

2° 

-.13131 

-.13020 

-.08904 

-.06304 

-.04638 

4° 

-.13249 

-.13540 

-.09223 

-.06775 

-.04964 

8° 

-.13738 

-.14780 

-.09900 

-.07536 

-.05761 

12° 

-.14305 

-.16186 

-.11036 

-.08804 

-.06812 

16° 

-.14932 

-.16660 

-.11813 

-.10181 

-.08188 

20° 

-.14266 

-.15690 

-.13008 



-.11775 

-.09565 

Table  A.9.  Tabular  data  for  CY 

1  Or 


a/M 

.8 

1.2 

2.0 

2.7 

3.5 

-8° 

.1635 

.1386 

.06877 

.04874 

.04270 

-4° 

.1639 

.1331 

.06996 

.05054 

.03745 

0° 

.1649 

.1357 

.07036 

.05269 

.04007 

2° 

.1608 

.1408 

.07194 

.05305 

.04232 

4° 

.1535 

.1474 

.07292 

.05556 

.04719 

8° 

.1452 

.1494 

.07273 

.06201 

.05618 

12° 

.1467 

.1553 

.07668 

.07168 

.06780 

16° 

.1414 

.1498 

.08360 

.08208 

.08052 

20° 

.1101 

.1231 

.09368 

.09570 

.08951 
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Table  A. 10.  Tabular  data  for  Cu 


a/M 

.8 

1.2 

2.0 

2.7 

3.5 

-8° 

.22653 

.35501 

.19609 

.16015 

.10417 

-4° 

.20408 

.26612 

.17228 

.12767 

.09510 

0° 

.13163 

.12245 

.08447 

.04391 

.02323 

2° 

.07857 

.03469 

.01228 

-.00701 

-.01779 

4° 

-.02106 

-.06000 

-.05953 

-.06421 

-.05336 

8° 

-.17469 

-.27837 

-.18083 

-.17934 

-.10744 

12° 

-.24204 

-.33265 

-.22214 

-.22066 

-.12232 

16° 

-.24408 

-.30857 

-.22214 

-.22030 

-.11416 

20° 

-.20857 

-.27714 

-.19200 

-.19077 

-.10417 

Table  A.  11.  Tabular  data  for  C, 


a/M 

.8 

1.2 

2.0 

2.7 

3.5 

-8° 

-.13341 

-.16085 

-.07361 

-.06506 

-.05441 

-4° 

-.14801 

-.17360 

-.09268 

-.08207 

-.07059 

0° 

-.15265 

-.17539 

-.10111 

-.09057 

-.07353 

2° 

-.15044 

-.17204 

-.10288 

-.08725 

-.06949 

4° 

-.14314 

-.16957 

-.10067 

-.07874 

-.05956 

8° 

-.13119 

-.15459 

-.09224 

-.06617 

-.04265 

12° 

-.12345 

-.13714 

-.08115 

-.06174 

-.03860 

16° 

-.11062 

-.12371 

-.07251 

-.05989 

-.04669 

20° 

-.10111 

-.11655 

-.06829 

-.05767 

-.04853 
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Table  A.  12.  Tabular  data  for  C, 


n^r. 


a/M 

.8 

1.2 

2.0 

2.7 

3.5 

-8° 

.22168 

.10270 

-.05374 

-.02637 

-.02044 

-4° 

.25664 

.16052 

-.02687 

.01319 

.00365 

0° 

.31726 

.22611 

.00330 

.05275 

.03285 

2° 

.33186 

.24984 

.02379 

.06300 

.04380 

4° 

.33496 

.26882 

.03877 

.07692 

.05547 

8° 

.32743 

.25070 

.07048 

.10549 

.07372 

12° 

.34115 

.19935 

.08899 

.13077- 

.09270 

16° 

.36991 

.20453 

.09956 

.12161 

.12190 

20° 

.35973  .16915 

.13084 

.15824 

.15547 

Table  A.  13.  Tabular  data  for  C\ 


6P 


a/M 

.8 

1.2 

2.0 

2.7 

3.5 

-8° 

-.033468 

-.022174 

.001062 

.002825 

-.004925 

-4° 

-.040358 

-.030174 

-.002743 

-.002082 

-.000373 

0° 

-.049933 

-.035826 

-.006726 

-.008327 

-.009403 

2° 

-.051007 

-.036957 

-.007965 

-.009814 

-.010896 

4° 

-.050828 

-.037217 

-.009204 

-.010409 

-.011045 

8° 

-.049664 

-.032957 

-.010929 

-.011004 

-.010000 

12° 

-.050380 

-.025130 

-.011770 

-.011004 

-.007313 

16° 

-.055302 

-.023391 

-.010265 

-.009665 

-.009701 

20° 

-.037673 

-.014957 

-.014159 

-.015762 

-.015821 
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Table  A.  14.  Tabular  data  for  C, 


a/M 

.8 

1.2 

2.0 

2.7 

3.5 

-8° 

61.3158 

4.0263 

-15.2343 

-9.6414 

-10.8909 

-4° 

74.2982 

40.0000 

-1.7188 

-5.3386 

-4.6446 

0° 

86.3158 

79.4748 

9.0625 

0.0000 

2.2422 

2° 

92.4561 

99.2560 

13.3594 

.1593 

5.2853 

4° 

97.6316 

117.4617 

16.3281 

.9562 

7.6877 

8° 

91.5789 

151.6849 

16.0156 

3.1076 

3.6036 

12° 

86.8421 

179.9562 

15.6250 

4.3825 

0.0000 

16° 

88.0702 

198.4245 

21.4063 

6.0558 

-1.2813 

20° 

93.9474 

294.0919 

25.9375 

7.8884 

3.0430 

Table  A. 15.  Tabular  data  for  C> 


a/M 

.8 

1.2 

2.0 

2.7 

3.5 

-8° 

-9.2259 

-7.7656 

-3.9717 

-3.1351 

-2.0118 

-4° 

-9.2417 

-7.6484 

-3.9717 

-3.1583 

-2.0197 

0° 

-9.2338 

-7.7656 

-3.9152 

-3.1390 

-2.0039 

2° 

-9.1945 

-7.9375 

-3.8799 

-3.0579 

-2.0039 

4° 

-9.1866 

-8.0937 

-3.8233 

-2.9730 

-2.0039 

8° 

-9.1945 

-8.3281 

-3.7385 

-2.8108 

-1.8619 

12° 

-9.2259 

-8.3203 

-3.6396 

-2.7027 

-1.7357 

16° 

-9.1081 

-8.2969 

-3.5830 

-2.6718 

-1.6331 

20° 

-9.0688 

-7.9609 

-3.6254 

-2.6873 

-1.6095 
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Table  A.  16.  Tabular  data  for  Cn 


ot/M 

.8 

1.2 

2.0 

2.7 

3.5 

-8° 

121.7839 

121.3270 

44.3658 

28.8446 

16.7480 

-4° 

98.7909 

94.1548 

41.3766 

27.5697 

16.7480 

0° 

76.4321 

66.1137 

37.2075 

26.7729 

16.4228 

2° 

64.9356 

51.5798 

34.2183 

25.9761 

15.9450 

4° 

53.1219 

37.9147 

30.9931 

25.0996 

15.7317 

8° 

30.2081 

9.3997 

24.0708 

22.3904 

13.4959 

12° 

7.6115 

-18.1675 

16.9912 

21.0359 

13.0081 

16° 

-15.7780 

-46.9984 

10.9341 

20.31-87 

12.8455 

20° 

-38.5332 

-77.0142 

6.3717 

20.1594 

12.7642 

Table  A.  17.  Tabular  data  for  C 


No 


M 

CNq 

.8 

-149.068 

1.2 

-134.226 

2.0 

-87.602 

2.7 

-77.176 

3.5 

-75.149 

Table  A.  18.  Tabular  data  for  C„ 


M 

Cmq 

.8 

-1981.268 

1.2 

-1900.576 

2.0 

-1218.300 

2.7 

-1040.346 

3.5 

-1007.925 
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Table  A.  19.  Tabular  data  for  C 


N* 


M 

CN. 

.8 

163.215 

1.2 

198.109 

2.0 

168.132 

2.7 

172.199 

3.5 

164.634 

Table  A.20.  Tabular  data  for  C„ 


M 

^m° 

.8 

442.575 

1.2 

532.153 

2.0 

94.406 

2.7 

15.694 

3.5 

-27.525 

Table  A.21.  Tabular  data  for  CYr 


M 

cYr 

.8 

219.964 

1.2 

186.303 

2.0 

92.659 

2.7 

53.447 

3.5 

40.197 
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Table  A. 22.  Tabular  data  for  Cnr 


M 

Cnr 

.8 

-1857.534 

1.2 

-1854.110 

2.0 

-1182.192 

2.7 

-937.671 

3.5 

-849.315 

Table  A.23.  Tabular  data  for  & 


M 

ClD 

.8 

-91.499 

1.2 

-114.634 

2.0 

-67.057 

2.7 

-47.486 

3.5 

-35.671 

APPENDIX  B 
INERTIAL  DATA  FOR  THE  EMRAAT  AIRFRAME 

P  A  It8  aPPfndi.^  §ives  the  weight  and  moments  and  products  of  inertia  of  the  EM- 
K.AA1  missile  with  no  fuel. 


Missile  Weight 
Moments  of  Inertia 


Products  of  Inertia 


W  =  227  lb 

Ixx  =  1.08  slug-ft2 
Iyy  =  70.13  slug-ft2 
Izz  =  70.66  slug-ft2 

Ixy  =  0.274  slug-ft2 

hz  =  0.704  slug-ft2 

Iyz  =  -0.017  slug-ft2 
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